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Abstract 

We introduce a method to reconstruct an element of a Hilbert space in terms of an arbi- 
trary finite collection of linearly independent reconstruction vectors, given a finite number of 
its samples with respect to any Riesz basis. As we establish, provided the dimension of the 
reconstruction space is chosen suitably in relation to the number of samples, this procedure 
can be numerically implemented in a stable manner. Moreover, the accuracy of the resulting 
approximation is completely determined by the choice of reconstruction basis, meaning that the 
reconstruction vectors can be tailored to the particular problem at hand. 

An important example of this approach is the accurate recovery of a piecewise analytic 
function from its first few Fourier coefficients. Whilst the standard Fourier projection suffers 
from the Gibbs phenomenon, by reconstructing in a piecewise polynomial basis, we obtain an 
approximation with root exponential accuracy in terms of the number of Fourier samples and 
exponential accuracy in terms of the degree of the reconstruction function. Numerical examples 
illustrate the advantage of this approach over other existing methods. 

1 Introduction 

Suppose that H is a separable Hilbert space with inner product (•, •) and corresponding norm ||-||. 
In this paper, we consider following problem: given the first m samples {{f,if>j)}JLi of an element 
/ 6 H with respect to some Riesz basis {V'jj^Li of H (the sampling basis), reconstruct / to high 
accuracy. Not only does such a problem lie at the heart of modern sampling theory [3TJ [37] , it also 
occurs in a myriad of applications, including image processing (in particular, Magnetic Resonance 
Imaging), and the numerical solution of hyperbolic partial differential equations (PDEs). 

In practice, straightforward reconstruction of / may be achieved via orthogonal projection with 
respect to the sampling basis. Indeed, for an arbitrary / € H, this is the best possible strategy. 
However, in many important circumstances, this approximation converges only slowly in m, when 
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measured in the norm on H, or not at all, if a stronger norm (for example, the uniform norm) is 
considered. 

A prominent instance of this problem is the recovery of a function / : T — > R from its first 
m Fourier coefficients (here T = [—1,1) is the unit torus). In this instance, H = L 2 (T) is the 
space of all 2-periodic, square-integrable functions of one variable. Provided / is analytic, it is 
well-known that its Fourier series (the orthogonal projection with respect to the Fourier basis) 
converges exponentially fast. However, whenever / has a jump discontinuity, its Fourier series suffers 
from the well-known Gibbs phenomenon [41] , Whilst convergence occurs in the L 2 norm, uniform 
convergence is lacking, and the approximation is polluted by characteristic O (1) oscillations near the 
discontinuity. Moreover, the rate of convergence is also slow: only O{ro~ 2 ) when measured in the L 2 
norm, and O (m _1 ) pointwise away from the discontinuity. Needless to say, the Gibbs phenomenon 
is a significant blight of many practical applications of Fourier series |35) . It is a testament to its 
importance that the design of effective techniques for its removal remains an active area of inquiry 

Returning to the general form of the problem, let us now suppose that some additional infor- 
mation is known about the function /. For example, / may be sparse in a particular basis (e.g. a 
polynomial of low degree) or may possess certain regularity. In particular, in the Fourier instance, 
we may know that / is piecewise analytic with jump discontinuities at known locations in T. In this 
circumstance, it seems plausible that a better approximation to / can be obtained by expanding 
in a different basis (e.g. a piecewise polynomial basis). To this end, let us introduce the so-called 
reconstruction space (of dimension n) and seek to approximate / by an element f n ,m consisting of 
n linearly independent elements of this space. 

As we will show in due course, provided reconstruction is carried out in a certain manner, a 
suitable approximation /„. m can always be found. Essential to this approach is that m (the number 
of samples) is chosen sufficiently large in comparison to n (equivalently, n is chosen sufficiently 
small in comparison to m). However, provided this is the case, the approximation f n>m inherits 
the principal features of the reconstruction space. In particular, /„ !?n is quasi-optimal (or, under 
certain conditions, asymptotically optimal), in sense that the error \\f — fn,m\\ can be bounded by 
a constant multiple of ||/ — Q n /||j where Q n f is the orthogonal projection onto the reconstruction 
space (in other words, the best approximation to / from this space). Moreover, from a practical 
standpoint, this method can be implemented by solving a linear least squares problem. Whenever 
the reconstruction vectors are suitably chosen (e.g. if they form a Riesz basis), the corresponding 
linear system is well-conditioned and the least squares problem can be solved in O (mn) operations 
by standard iterative techniques. 

Consider once more the example of Fourier series. Let / : T — ¥ R be an analytic function with 
jump discontinuity at x = — 1 (equivalently, an analytic, nonperiodic function). As mentioned, the 
Fourier series of / lacks uniform convergence. However, since / is analytic, it makes sense to seek 
to reconstruct / in a system of polynomials. It is well-known that the best n th degree polynomial 
approximation of an analytic function converges exponentially fast in n As we shall prove, 
with n = O {y/m), we obtain a quasi-optimal polynomial approximation /„ jm to / from only its 
first m Fourier coefficients, regardless of the particular family of polynomials used. This results in 
exponential convergence of f n ,m to / in n (the polynomial degree), or root exponential convergence 
in m (the number of Fourier samples). Moreover, whenever Legendre polynomials are used, the 
approximation is computable in a stable manner in only O (mn) operations. The use of, for example, 
Chebyshev polynomials results in a method with O (n) condition number, that can be implemented 
in O(mnJ) operations. Furthermore, whilst retaining the aforementioned features, this procedure 
can be easily generalised to recovery of a piecewise analytic function of one variable (using piecewise 
polynomial bases), and to the case of multivariate functions defined in tensor-product regions. 

There are a number of existing algorithms for the removal of the Gibbs phenomenon from Fourier 
series. One of the most well-known, which also provides a polynomial approximant, is the Gegen- 
bauer reconstruction technique 133 EH] ■ As we discuss further in Section|3| the method developed 
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in this paper has a number of key advantages over this device. Numerical results also indicate its 
superior performance for a number of test problems. 

The method developed in this paper was previously introduced by the authors in [3] within the 
context of abstract sampling theory. Whilst this problem, in this abstract form, has been extensively 
studied in the last couple of decades (in particular, by Eldar et al [201121]) see also 07]), to the best 
of our knowledge this particular method does not appear in any existing literature. For a more 
detailed discussion of the relation of this approach to existing schemes we refer the reader to [3]. 
Conversely, in this paper, after presenting the general version of the method in abstract terms, 
we will focus primarily on its application to the Fourier coefficient reconstruction problem. On 
this topic, a similar approach, but only dealing with reconstructions in Legendre polynomials from 
Fourier samples of analytic functions, was discussed in [32 . This can be viewed as a special case of 
our general framework. Furthermore, by examining this example as part of this framework, we are 
able to extend and improve the work of [32] in the following ways: (i) we derive a procedure allowing 
for reconstructions in any polynomial basis, not just Legendre polynomials, (ii) we extend this 
approach to reconstructions of piecewise smooth functions using (arbitrary) piecewise polynomial 
bases, (iii) we generalise this work to smooth functions of arbitrary numbers of variables and (iv) 
we obtain improved estimates for both the error and the necessary scaling n = O (s/m) required for 
implementation. 

Aside from yielding these improvements, a great benefit of the general framework presented in 
this paper is that it is immediately applicable to a whole host of other reconstruction problems. To 
illustrate this generality, in the final part of this paper we consider its application to the accurate 
reconstruction of a piecewise analytic function from its orthogonal polynomial expansion coefficients. 
Such a problem is typical of that occurring in the application of polynomial spectral methods to 
hyperbolic PDEs [27], where the shock formation inhibits fast convergence of the polynomial approx- 
imation. As we highlight, this issue can be overcome in a completely stable fashion by reconstructing 
in a piecewise polynomial basis. 

The outline of the remainder of this paper is as follows. In Section[2]we introduce the reconstruc- 
tion procedure and establish both stability and error estimates. Section [3] is devoted to (piecewise) 
polynomial reconstructions from Fourier samples. In Section [4] we consider reconstructions from 
tensor-product spaces, and in Section [5] we discuss other reconstruction problems. Finally, in Sec- 
tion [6] we present open problems and challenges. 



2 General theory of reconstruction 



In this section, we describe the reconstruction procedure in its full generality. To this end, suppose 
that {V'jljli is a Riesz basis (the sampling basis) for a separable Hilbert space H over the field 
C. Let (•,•) be the inner product on H, with associated norm ||-||. Recall that, by definition, 
span{i/>i, ip2, ■ ■ ■} is dense in H and 



3=1 



3=1 



*3l > 



Va = {ai,a 2 ,...} € / 2 (N), 



(2.1) 



3=1 



for positive constants ci,C2. Equivalently, ipj = U(^j), where {^jl^Li is an orthonormal basis for 
H and U : H — » H is a bounded, bijective operator. Using this definition, it is easy to deduce that 
{^}^i also satisfies the frame property 



di||/|| 2 <EK/'^>l 2 ^2||/|| 2 , V/6H, 

3 = 1 



(2.2) 



for d\,di > 0, where the largest possible value for d,2 is H^/Hf^jj and the smallest possible value for 

di is Hw-Ih^h D3|- 
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Suppose now that the first m coefficients of an element / € H with respect to the sampling basis 
are given: 

fj = (f,i>j), J=l m. (2.3) 

Set S m = spanj^i, . . . , i/j m } and let V m : H — > S m be the mapping 

m 

f^V m f = J2(f>^)^r (2-4) 

3=1 

We now seek to reconstruct / in a different basis. To this end, suppose that {<j>i, . . . , cf> n } are linearly 
independent reconstruction vectors and define T„ = span{</>!, . . . , <p n }. Let Q n : H — > T„ be the 
orthogonal projection onto T„. Direct computation of Q n f, the best approximation to / from T n , 



is not possible, since the coefficients (f,<fij) are unknown. Instead, we seek to use the values (2.3) 
to compute an approximation f n . m € T„ that is quasi- optimal, i.e. ||/ — Q n f\\ < 11/ — fn,m\\ < 
C\\f — Qnf\\ for some constant C > independent of /, n and m. To do this, we introduce the 
sesquilinear form a m :HxH->C, given by 

a m (g, h) = (V m g, h) , V<?, h e H. (2.5) 

Note that, since 

rn 

(P m g, h) =Y^ (f, ^j) (g, V'j) = a m{g, /), V/, g e H, 

3=1 

a m is a Hermitian form on H x H (here z is the complex conjugate of z € C). With this to hand, we 
now define /„ im by the following condition 

a m (fn, m , 4>) = a m (f, 4>), V0 e T„. (2.6) 

Upon setting = cj)j, j = l,...,n, this becomes an n x n linear system of equations for the 
coefficients ai, . . . , a n of f n ,m = Y^i=i a j ( t'j- We shall defer a discussion of the computation of this 



approximation to Section 2.3 in the remainder of this section we consider the analysis of /„ 



Before proving the main theorem regarding (2.6), let us first give an explanation as to why this 



approach works. As mentioned, key to this technique is that the parameter m is sufficiently large in 



comparison to n. To this end, let n be fixed and suppose that m — » oo. Due to (2.1 ), the mappings 
V m converge strongly to a bounded, linear operator V , given by 

oo 

Vf = J2(f^i)^> V/eH. (2.7) 

3=1 



Hence, for large m, the equations (2.6) defining /„ m resemble the equations 



a(/ n , <£) = «(/, V^T„, (2.8) 

where a : H x H — > C is the Hermitian form a(f,g) = (Vf,g)- Thus, it is reasonable to expect that 
fn,m ~^ fn as Tn — Y oo , provided such a function f n exists. However, 

Theorem 2.1. For all n£N, the function f n exists and is unique. Moreover, 

||/-/n||<|||/-Q„/||, (2.9) 



where di and di arise from (2.2) 
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This theorem can be established with a straightforward application of the Lax-Milgram theorem 



and its counterpart, Cea's lemma [25]. Indeed, due to (2.2) and (2.7), 

oo 

di\\g\\ 2 < a(g,g) = £ | (g,^) | 2 < d 2 ||g|| 2 , Wg e H. 
i=i 



(2.10) 



Hence the form a(-, ■) defines an equivalent inner product on H. Nonetheless, we shall present a 
self-contained proof, since similar techniques will be used subsequently. 

Proof. Let U : T„ — > C n be the linear mapping g i-> {(Vg, 4>j)}$=v To prove existence and unique- 
ness of f n it suffices to show that U is invertible, upon which it follows that /„ = W _1 {("P/, </>j)}™ =1 . 
Suppose that lAg — 0. Then, by definiti on, (V g, 4>j) — for j = 1, . . . , n. Using linearity, we deduce 
that (Vg,g) = 0. Now, it follows from ( |2.10[ ) that = (Pg,g) > di||.a|| 2 , giving 5 = 0. Hence, U is 
invertible and /„ exists and is unique. 

Now consider the error estimate (2.9). Using ( |2.2[ ) once more, we obtain 

1 



ii/ -Mr <j-(vu-u),f-u 



By definition of /„, (P(f - f n ),4>) = 0, V<j> € T n . In particular, setting </> = /„- Q n /, yields 



II/- /nil 2 < ~[- \P(f- fn),f~ Qnf 

Since a(-, •) forms an equivalent inner product on H, an application of the Cauchy-Schwarz inequality 
gives 

11/ - /„|| 2 < j~ \{V(f - /„), / - f n )(V(f - Q n f), f - Q n /)1 * < ^-\\f - /„||||/ - Q n f\\, 



di 



di 



as required. 



□ 



This theorem establishes existence and quasi-optimality of /„ « f n ,mi thereby giving an intuitive 
argument for the success of this method. We now wish to fully confirm this observation. To this 
end, let 



C n>m = inf V" | (<f>,ipj) | 2 = inf (V m (j), 4>) = inf a m (</), i 

06T„ 06T„ 06T„ 

11011=1 J=1 11011 = 1 11011 = 1 

The quantity C„. m plays a fundamental role in this paper. Note that 



(2.11) 



Lemma 2.2. For all n,meN,0< C„. m < cfe- Moreover, for each n, C„,m — > C* > d\ as m — > oo, 
where 

oo 

06 Ire 06 In 06 In 

11011 = 1 3 11011 = 1 11011 = 1 



and c?i is defined in (2.2). 

Proof. Consider first the quantity 



£n,m = SUp I I 2 = SU P (P<l > -'Pm<f>, 



*6Tn 



06T„ 



(2.12) 
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Due to (2.2), the infinite sum is finite for any fixed <fi, and tends to zero as m — > oo. Now let 

be an orthonormal basis for T„ and set <f> — Ej=i a 3^3- Two applications of the Cauchy-Schwarz 

inequality gives 

n 

Eito^> i a <iwi a EEi<**>^>i 2 - 

j>m k—lj>m 

Hence e n , m < X)fc=i X^> m I (^kiipj) | 2 j an d we deduce that e niTO is both finite and e n>m — > as 
to — > oo. Noticing that |C„. m — C*| < e„ im , Vn, m G N, gives the first part of the proof. For the 
second, we merely use (2.2 1. □ 



Aside from C n>m , we also define the quantity D n>m by 



sup sup |(P m /,5)| . 

i|/||=l llflll=l 



(2.13) 



For this, we have the following lemma: 



Lemma 2.3. For all m, n 6 N, < F) n ^ m < d 2 - Moreover, suppose that V is such that V (T„) C T„ 
(7or example, when V = I is the identity), then D„ m < C2e„ im , where e n , m is as in (2.12). In 
particular, for fixed n, D n m — » as to — > oo. 

Proof. Let /, <? € H. By definition, 



3=1 



3=1 



E i V'j 

3=1 



(2.14) 



Hence (2.2) now gives the first result. Now suppose that / g T„ and V(T n ) C T„. Since P m is 
self-adjoint, we have (V m f,g) — (fiV m g) — (f,V m g ~ Vg)- Here the second equality is due to the 
fact that / _L T„ and Vg G T„ for g G T„. By the Cauchy-Schwarz inequality, we obtain 

D n , m < sup \\Vg-V m g\\. 

<?6T„ 

II a ll = i 



For g G T„, we note from (2.1) that 



\Vg - V m gf <c 2 J2\ (g, I 2 = c 2 (Vg - V m g, g) < c 2 e r , 



where the final equality follows from the definition of e„.. 



□ 



D„ 



and C* 



At this moment, we mention the following important fact. The constants C n ^ m , j^ n ,m 
as well as the approximation f nm , are determined only by the space T„, not by the choice of 
reconstructions vectors 4>i, - ■ ■ ,4>n themselves. As we shall discuss later, the choice of reconstruction 
vectors only affects the stability of the scheme. 

This observation aside, we are now able to state the main theorem of this section: 

Theorem 2.4. For every n G N there exists an too such that the approximation f n ,m, defined by 
\2. 6V, exists and is unique for all to > mo, and satisfies the stability estimate \\f n ,m\\ < d 2 C~ ^11/11- 
Furthermore, 



\\f-Qnf\\ < ||/-/n,m|| < K n , m \\f - Q n f\\, K n , m = y/ 1 + D^C^m 

Specifically, the parameter Too is the least value of to such that C nm > 0. 



(2.15) 
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To prove this theorem, we first recall that a Hermitian form a : H x H — ► K is said to be continuous 
if, for some constant 7 > 0, \a(f, g)\ < 7||/|| ||.g| for all f,g £ H. Moreover, a is coercive, provided 
d(f, /) > w II/I| 2 j V/ £ H, for w > constant [25]. We now require the following lemma: 

Lemma 2.5. Suppose that a m : H x H — > K is i/ie sesquilinear form a m (f,g) — (V m f,g). Then, 
a m is continuous with constant 7 < di . Moreover, for every n £ N there exists an mo such that the 
restriction of a m to T„ x T„ is coercive for all m > mo- Specifically, if C„ jm is given by (2.11), 
then mo is the least value of m such that C n , m > 0, and, for all m > mo, a m (f,f ) > C n . m \\f\\ 2 , 
V/ £ T„. Finally, for all f £ H and g £ T„, we have a m (f - Q n f,g) < D n , m \\f - Qn/HIMI- 



the definitions (2.111 and (2.13) of C. 



Proof. Continuity follows immediately from (|2.14|). For the second and final results, we merely use 

respectively. 



□ 



Proof of Theorem \2.4\ To establish existence and uniqueness, it suffices to prove that the linear 
operator U : T„ — > C", g H> {(V m g, ^j)}" =1 is invertible. Suppose that g £ T n with Ug = 0. By 
definiti on, w e have (V m g,<j>j) = for j = 1, . . . , n. Using linearity, it follows that (V m g,g) = 0. 
Lemma 2.5 now gives < C nim ||g|| 2 < 0. Hence g = 0, and therefore U has trivial kernel. 

Stability of f n ^ m is easily established from the continuity and coercivity conditions. Setting 
<t> = fn,m in J2l3|) gives 



Cfi,m 1 1 fn,m || — ^rn {fn.m ? fn.rn ) — &m (</"; fn,m) — ^2 \ \f\ \ \ \ fn,m \ \ 



a m (f - Qnf, 



as required. Now consider the error estimate (|2.15|). Suppose that we define e n _ rn — f n , m ~Qnf £ T n . 
Then, by definition of /„, m , we have a m (e n ^ 
4> = e n , m , we obtain 

II ^n,m | 

en,m|| 2 ~T~ || J — J 

□ 



|/-Qn/||. 

Since Q„f is the orthogonal projection onto T„, we have ||/ — / n ,m|| 2 
which gives the full result. 



£ T„. In particular, setting 
(2.16) 

1/ - Qn/H 2 , 



When a m is shown to be continuous and coercive, it may be tempting to seek to apply the 



Lax-Milgram theorem and Cea's lemma to obtain Theorem 2.4 However, the Hermitian form a m , 
when considered as a mapping H x H 4 C, will not, in general, be coercive. This is readily seen 
from the definition of C„ jm . The finite-dimensional operator V m \T n converges uniformly to V\t„, 
whereas its infinite-dimensional counterpart V m typically does not (for example, when V m is the 
Fourier projection operator and H = L 2 (T)). Hence, a m only becomes coercive when restricted to 
T„ x T„, and these standard results do not automatically apply. 



Though Theorem 2.4 establishes an estimate for the error / — / n , m measured in the natural norm 



on H, it is also useful to derive a result valid for any other norm defined on a suitable subspace of 
H (for example, this may be the uniform norm on T in the case of Fourier series). To this end, let 
I • (I be such a norm and define G = {g £ H : \\g\\ < oo}. We have 



Corollary 2.6. Suppose that f £ G, T„ C G and that /„, 

771 > 7770, 

1/ - fn,ml < If ~ Gn/| + ^^||/ - Qn/I 



n 



is defined by (2.6). Then, for all 

(2.17) 



where k„ 



sup^gx,, 

11011=1 



and C. 



a. m> T) nm 



Proof. Let e n>m = fn,r. 



are given by (2.11) and (2.13) respectively. 
Q n f once more. Since e nm £ T„, it follows from the definition of k n and 



the inequality (2.16) that 



k D 

|pn,m||| — k n \\e njm 

II < -^^Wf-QnfW 



The full result is obtained from the triangle inequality 



fn,r 



< 



- Qn/I- 



□ 
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Remark 1 In practice, it is useful to have an upper bound for the constant k n . A simple exercise 
gives k n < Y^j=i ll^illi where {$j}" = i is any orthonormal basis for T n . 



Theorem 2.4 confirms that the approximation f n . m is quasi-optimal whenever m is sufficiently 
large in comparison to n. In other words, the error ||/ — f n .m\\ can be bounded by a constant 
multiple (independent of n) of ||/ — Q n f\\. Naturally, whenever {</> l7 . . . , <f> n } are the first n vectors 
in an infinite sequence <pi , fa , ■ • ■ that form a basis for H (a natural assumption to make from a 
practical standpoint), we find that f n ,m converges to / at the same rate as Q n f, provided m scales 
appropriately with n. Under the same assumptions, Corollary |2.6| also verifies convergence of f n ,m 
to / in I • II, whenever Q n f — > / in this norm and k n \\f — Q n f\\ — > as n — > oo. 

Note on other feature of this framework. Provided m is such that C nm > 0, any function / £ T„ 
will be recovered exactly by f n . m . In the language of sampling theory, this property is commonly 
referred to as perfect reconstruction. 

Remark 2 Whilst this framework relies on the fact that m can range independently of n, a natural 
question to ask is what happens if m is set equal to n (note that, in this case, one recovers the 
well-known consistent sampling framework of Eldar et al |21l 147] ). This question was discussed in 
detail in [3J, where it demonstrated that such an approach often leads to severe ill-conditioning as 
n = m — > oo. Additionally, stringent restrictions are placed on the types of vectors / that can 



be reconstructed (see also Section 3.6). Conversely, by allowing m to vary independently of n, we 
obtain a reconstruction f nm that is guaranteed to converge for any vector / £ H. Moreover, as we 
discuss in Section |2.3| provided the reconstruction vectors are suitably chosen, the computation of 
f n . m is completely stable. 

2.1 The case of an orthogonal sampling basis 

In the previous setup, the sampling basis was assumed to be a Riesz basis. The particular case 
of {tfij^^-i being an orthonormal basis warrants further study, since this situation often arises in 
applications (e.g. Fourier sampling). 



In this setting, both (2.1) and (2.2) (which are now identical) hold with equality and with all 



constants being precisely one. In other words, Parseval's identity 

oo 

ii/ii 2 =Ek/'^>i 2 > v / £H ' 

holds. Our first result provides a geometric interpretation for the constant C„ im in terms of subspace 
angles. Recall that if U and V are subspaces of H, then the angle #uv is given by 

cos(# uv )= inf ||Cvt*||, (2-18) 

IMI=i 

where Qv : H — > V is the orthogonal projection. We have 

Lemma 2.7. Suppose that {"0j}^i * s an orthonormal basis of H. Then C n . m — cos 2 (#), where 
8 = 0T n s m is the angle between the subspaces T„ and S m . 

Proof. Whenever {tpj}fLi is orthonormal, the operator V m is the orthogonal projection onto S m . 



Hence the result follows immediately from (2.11) and (2.18). □ 



Whilst this property is interesting, the following result has much greater significance: 

Theorem 2.8. Suppose that is an orthonormal basis of H. Then, for all m > mo, the 

approximation f n>m satisfies 

11/ - Qnf\\ < ||/ - fn,m\\ < K n , m \\f - Q n f\\ (2.19) 
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where 

K n , m = \Jl + (1 - C„, m )GT,m = \/l + tan 2 6>sec 2 6>, 



and 9 is as in Lemma 2.7 In particular, for fixed n, K nm — > 1 asm-) oo, and hence f nm — > Q n f 
as m —¥ oo. 



Proof. Note that V =1 for an orthonormal sampling basis. In particular, e njm = C* — C„ >m and 
C* = 1. Hence, using Lemma 2.3 we find that D 2 m < 1 — C n . m . The result now follows from 
Theorem [231 □ 



From this we conclude the following: not only is f n ,m quasi-optimal, it is also asymptotically 
optimal in the sense that f n , m —> Qnf, the best approximation to / from T„, asm-} oo. Hence, 
using this approach, we can recover an approximation to / that is arbitrarily close to the error 
minimising approximant (which, as mentioned, cannot be computed directly from the given samples). 
Moreover, the rate of convergence of f n ,m to Q n f is completely independent of the particular vector 
/, and relies only on rate of decay of the parameter 1 — C„, m . 

Note that asymptotic optimality also occurs for general Riesz bases whenever V(T n ) C T„, in 



which case D n m — ¥ as m — > oo (see Lemma 2.3 1. The case of orthonormal sampling vectors 



presents the most obvious example of a basis satisfying this condition. 

Remark 3 Whenever the vectors {4>j}JL\ are not orthonormal, a natural question to ask is whether 
we can modify the approach for computing / n , m to recover asymptotic optimality. This can be easily 



done, at least in theory, by replacing the operator T> m , given by (2.4), with the orthogonal projection 



H — > T„. In this case, both Lemma 2.7 and Theorem 2.8 will hold for nonorthogonal sampling bases. 
The downside of the approach is that it requires additional computational cost to compute f n ,m, as 
we explain at the end of the next section. 

Another potential means to recover asymptotic optimality is to define V m g = Ylj—i (g,4'j) ipj, 
where {ipj} is the set of dual vectors to the sampling vectors {ipj}. In this case, V m — > I strongly, 
and asymptotic optimality follows. In practice, however, one may not have access to the dual vectors, 
thus this approach cannot necessarily be readily implemented. 

2.2 Oblique asymptotic optimality 

Whenever the sampling basis is orthonormal, f n ,m converges to the best approximation Q n f as 
m — > oo. An obvious question to ask is what can be said in the general case, where the vectors {'ipj} 
only form a Riesz basis? The intuitive explanation given previously indicated that f nm ~ /„, where 



/„ is defined by (2.8). We now wish to confirm this observation. In fact, as in the orthonormal case, 
we may demonstrate a stronger result: namely, for fixed n € N, f n .m — > fn as m — > oo at a rate that 
is independent of the particular vector / G H. 

Recall that the form a(-, •) yields an equivalent inner product on H. Since /„ is defined by the 
equations a(f n ,<p) — a(f,<f>), V0 £ T n , the mapping / i-> /„ is the orthogonal projection onto T n 
with respect to this inner product. Letting \\g\\ a — y/a(g, g) be the corresponding norm on H, we 
now define the constants 

C n . m = inf {V m (f>,4>), Ai,m= sup sup \{V m f,g)\. (2.20) 

0eT„ , t _l geT 

IWL=i ||/f^illffl|.=i 

In this instance, T^ is defined with respect to the a-inner product, i.e. T^ = {/ € H : a(f,4>) = 
0, € T„}. Conversely, when considered with respect to the canonical inner product, this subspace 
is precisely P(T„)^ = {/ eft : (f,<f>) = 0, G V(T n )}. 



Note the similarity between C n ^ m and D n ,m and the quantities C n ,m and D n ^ m defined in (2.11 ) 



and (2.13) respectively. Roughly speaking, the former measure the deviation of f nm from Qnf, 



whereas, as we will subsequently show, the latter determine the deviation of f nm from /„. 
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With these definitions to hand, identical arguments to those given in the proofs of Lemmas 2.2 



and 2.3 now yield: 



Lemma 2.9. For all is as in X2.11V. Moreover, for fixed n, 
C n , m — > 1 as m — > oo. 

Lemma 2.10. For all m,neN, D n m < di and £>„ m < 02(1 — C n ,m)- I n particular, for fixed n, 
D n ,m ~~ as m — > oo. 

Using these lemmas, we deduce 



Corollary 2.11. If f n>m and f n are given by (2.6) and (2.8) respectively, then 



||/n,m /nlla — ~ 11/ /n||ai 



and we have the error estimate 



11/ — /n|| a < 11/ — /n,m||a < ^n,m||/ — fn\\a, K-n,m — \j 1 + ^\\,mPn%i- 

In particular, for any f £ H ; fn.m — ^ fn — y oo. 
Proof. Since (/„, m - /„) € T„, we have 

Cn,m 1 1 fn.m /n||a — (^m^fn.rn fn ) ; fn,m fn\ • 

Moreover, because (P m fn,m,<f>) — {V m f,4>), we deduce that 

Cn , m 1 1 fn . m fn 1 1 a — ^ !~*m ( / fn ) ; fn . m fn \ — Dn , m 1 1 / fn 1 1 a | /n , m /n 1 1 a j 



where the second inequality follows from the definition ( 2.20 1 of Z) n ,m and the fact that (f—f n ) <= 



the orthogonal complement of T„ with respect to the a-inner product. □ 

Note that the mapping W„ : / i— > /„ is an oblique projection with respect to the inner product 
(•, •) on H. In particular, W„ has range T„ and kernel ■p(T n )- L , and we have the decomposition 
H = T„ ©P(T„) ± . For this reason, we say that f n ,m possesses oblique asymptotic optimality. 

2.3 Computation of / n m 

Recall that the computation of the approximation /„, m involves solving the system of equations 



(2.6). This can be interpreted as the normal equations of a least squares problem. Suppose that 



fn,m = Z)i=i a i4>ji a = ("l, • ■ ■ ! d n ) £ C" and / = (/i, . . . , f m ). If U is the m x n matrix with 



(j,k) th entry (<j>k,ipj), then (2.6) is given exactly by Aa = Wf, where A = U^U and W is the 
adjoint of U. In other words, the vector a is the least squares solution of the problem Ua rj /. 

This system can be solved iteratively by applying conjugate gradient iterations to the normal 
equations, for example. The number of required iterations is dependent on the condition number 
k(A) of the matrix A. Specifically, the number of iterations required to obtain numerical convergence 
(i.e. to within a prescribed tolerance) is proportional to \J k(A) [26]. In particular, if k(A) is O (1) 
for all n and m > mo, then the number of iterations is also O (1) for all n. Hence, the cost of 
computing f n . m is determined solely by the number of operations required to perform matrix-vector 
multiplications involving U . In other words, only O (mn) operations. 
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Naturally, aside from this consideration, the condition number of A is also important since it 
determines susceptibility of the numerical computation to both round-off error and noise. Specifi- 
cally, an error of magnitude e in the inputs (i.e. the samples fj,j = l,..., m) will yield an error of 
magnitude roughly K,(A)e in the output f n , m - 

For these reasons it is of utmost importance to study the condition number of A. For this, we 
first introduce the Hermitian matrix A € c™ x ™ with (j,k) th entry (cf)j,<f>k). Note that A is the 
Gram matrix of the vectors {(f>%, ■ ■ ■ , </vi}- In particular, k(A) is a measure of the suitability of the 
particular vectors in which to compute Q n .f ■ With A to hand, we also introduce the related matrix 
A a € C" x ™ with (j, k) th entry a(0j, (f>k) — (V<frj , <j>k) , i- e - the Gram matrix with respect to the inner 
product a(-, •). 

The following lemma comes as no surprise: 

Lemma 2.12. The matrices A and A a are spectrally equivalent. In particular, for all n € N, 

^k(A) < n{A a ) < ^k(A). 

Proof. For any Hermitian matrix B, the condition number is the ratio of the largest and smallest 
eigenvalues (in absolute value). Moreover, if B is positive definite, then 

. „ [a^Ba^ , , _ . (a'Ba] ,_. 

inf <^ —— } = A min (B , sup i —— \ = A max (B . 2.21 
If (f> — X)i=i a j4>ji then a^ Aa = \\(f>\\ 2 and o^A a O! = a(</>, 0). Hence, spectral equivalence now follows 



immediately from (2.10 1. □ 



Concerning the condition number of the matrix A, we now have the following: 



Lemma 2.13. Suppose that m > m , where mo is as in Theorem 2.4, and C n , 



by (2.11) and (2.20) respectively. Then 



C n . m K{A a ) < K (A) < J—K(A a ), < k(A) < ^-«(i). 

Moreover, for fixed n, A — > A a as m — > oo, and, ifV=X, A — > A = A a . 
Proof. The m atrix A is Hermitian and, provided m > mo, positive definite. Hence, its eigenvalues 



are given by (2.21). For <j> — Y^j=i a j < i ) ii we have cv Aa = (V m 4>,<l>)- By definition of C n . m and 



C n , m , we find that A min (A) > C n . m X min (A) and X m in(A) > C' n , m \ min {A a ). Moreover, by (2.1) we 



have A max (A) < c?2A max (^4) and A max (j4) < A max (A ). The first result now follows immediately from 



(2.21 ). For the second, we merely note that each entry of A converges to the corresponding entry of 
A a as m — ¥ oo. □ 



Note the important conclusion of this lemma: computing f n ^ m from ( 2.6 ) is no more ill-conditioned 
than the computation of the orthogonal projection Q n f or the oblique projection W n f in terms of 
the vectors {(f>i, ...,</>„}. In practice, it is often true that these vectors correspond to the first n 
vectors in a basis {4>j}f^i of H with additional structure. Whenever this is the case, as the following 
trivial corollary indicates, we can expect good conditioning: 

Corollary 2.14. Suppose that {cf>j}'jL l is a Riesz basis for H (with respect to (•, •) ) with constants 
and c' 2 . Then 

-w < 



n 



Proof. This follows immediately follows from (2.1) and Lemma 2.13 □ 



Put together, a fundamental conclusion of Theorem 2.4 Lemma 2.13 and Corollary 2.14 is the 
following: for a given reconstruction space T n , the individual vectors </>i,...,0„ can be chosen 
arbitrarily, without altering the analysis of the approximation /„, m (which itself does not depend 
on the individual vectors used to represent it). The choice of vectors only becomes important when 
considering the condition number of linear system to solve. Moreover, the quality of a system of 
vectors for the reconstruction problem is completely intrinsic, in that it is determined only by the 
corresponding Gram matrix (in particular, it is independent of the sampling space). 



Corollary 2.14 confirms that the approximation f n ^ m can be readily computed in a stable manner 
for many choices of reconstruction basis. However, to fully implement this method, as we discuss 
further in the next section, it is useful to have numerical way of computing C n ^ m . The following 
lemma provides such a means: 

Lemma 2.15. The quantity C rhm is given by C n . m — \ min (A~ 1 A). Moreover, if A and A commute, 
then C n m = 1 — ||7 — .A _1 .A||. In particular, if {4>j}™ = i is an orthonormal basis, then C rlyTn — 
X mixi (A)'=l-\\I-A\\. 

Proof. By definition C nm = inf ^Pm4>i 4>). Letting cf> = X),-=i a j4>ji we fmd that 

11011 = 1 

„ . , Y,lk=i a j~ s k('P m (t>j,<t>k) . , otAa 

C n .m = mf _ , 7 7 \ = mf „ ■ 

We now claim that, for arbitrary Hermitian positive definite matrices B and C with B nonsingular, 
the following holds: 

mf t = X min (B C), sup —r— = X max {B C). 
aeC" a<Ba Q6C ™ a'Ba 

To do so, write B = D^D, with D nonsingular. Then, after rearranging, we obtain 

lnf +D = mf 7T75 = X min {D l CD ), 

for example. However, a trivial calculation confirms that the eigenvalues of D~ 1( CD~ 1 are identical 
to those of B~ 1 C, thus establishing the claim. Since A is nonsingular, this confirms that C nm — 
A m in(^4 _1 ^4)- For the second result, we merely notice that A m i n (£>) = 1 — A max (7 — B) = 1 — — B\\, 
whenever B is Hermitian. □ 

In Section [2.1 1 we briefly discussed a modified approach where the operator V m , usually given by 



(2.4), was replaced by the orthogonal projection operator. The advantage of this approach is that 
it guarantees asymptotic optimality. However, the downside is additional computational expense. 
Indeed, the corresponding matrix is of the form A = WV^ 1 !/, where V G C mxm has (j,k) th 
entry (i/}j,i/)k)- Hence, if conjugate gradients iterations are used, at each stage we are required to 
compute matrix-vector products involving the m x m matrix V^ 1 (assuming that V" 1 had been 
precomputed) . In general, this requires O (m 2 ) operations. Thus, we incur a cost of O (m 2 ), as 
opposed to O (mn) for the original algorithm. Hence, in practice it may be better settle for only 
quasi- and oblique asymptotic optimality, whilst retaining a lower computational cost. 
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2.4 Conditions for guaranteed, quasi-optimal recovery 



Let us return to the standard form of the method once more. To implement this method, it is 
necessary to have conditions that guarantee nonsingularity, stability and quasi-optimal recovery. In 
other words, for given sampling and reconstruction bases, we wish to study the quantity 

0(n; 6) = min {m G N : C n . m > 9} , 0e(O,d 2 ), (2.22) 



where C n , m is given by (2.111 and dn stems from (2.2| . Note that C n , m < &i by (2.2), thereby 



explaining the stated range of 8. Also, by Lemma 2.2 we have that lim^j^oo C n ^ m > di > 0, thus 
is well-defined. 

By definition, 6(n; 9) is the least m such that ||/ — f n ,m\\ < c (0)ll/ — Sn/||i where 

c{6) = yjl + dp- 2 or y/l + (1 - 0)0-2, (2.23) 

whenever the sampling basis is orthonormal. In other words, the least m required for quasi-optimal 
recovery with constant c(0). Thus, provided m > 0(n;0), the approximation /„ ilTl converges at 
the same rate as Q n f as n —> oo. In addition, m > 0(n;0) guarantees that ||/ n .m|| < rf2^ _1 ||/|| 
and n(A) < d20~ 1 n(A), thus making the linear system for f niTn solvable in a number of operations 



proportional to y d26~ x K{A). 

Note that 0(n; 0) is determined only by the sampling and reconstruction spaces S m and T n . 
Whilst 0(n;0) can be numerically computed for any pair of spaces via the expression given in 
Lemma |2.15| analytical bounds must be determined on a case-by-case basis. In the next section, 
where we consider the recovery of functions from their Fourier samples using (piecewise) polynomial 
bases, we are able to derive explicit forms for such bounds. 

Remark 4 As mentioned, the framework developed in this section was first introduced by the 



authors in [3j. Whilst a result similar to Theorem 2.4 was proved, there are a number of important 



improvements offered by the theory presented in this paper: 

1. In [3] it was assumed that the reconstruction vectors 4>i, . . . , <p n were the first n in an infinite 
sequence of vectors that formed a Riesz basis for H. Conversely, Theorem |2.4| depends only on 
the subspace T n , and thus the individual reconstruction vectors can be chosen arbitrarily. 

2. The constants K n m and C n , m are known exactly in terms of the sampling and reconstruction 
bases, and, whenever the sampling vectors are orthogonal, they possess a simple geometric 
interpretation in terms of subspace angles. Moreover, these constants can be computed by 
determining either the minimal eigenvalue or, in certain cases, the norm of an n x n matrix. 

3. Simple, explicit bounds for the condition number of the matrix A are known in terms of the 
constant C n ^ m and the Gram matrix A. 

4. The behaviour of f n , m as m — ^ oo (for n fixed) can be fully explained in terms of oblique 
asymptotic optimality. 



3 Polynomial reconstructions from Fourier samples 

One of the most important examples of this procedure is the reconstruction of an analytic, but 
nonperiodic function / to high accuracy from its Fourier coefficients. Direct expansion in Fourier 
series converges only slowly in the L 2 (— 1, 1) norm, and suffers from the Gibbs phenomenon near the 
domain boundary. Hence, given the first m Fourier coefficients of /, we now seek to reconstruct / 
to high accuracy in another basis. 

Let H = L 2 (— 1, 1) (the space of square integrable functions on (—1, 1)), / : (—1, 1) — > K and 
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be the standard Fourier basis functions. For m > 2, we assume that the coefficients 



fj = / f(x)ipj(x)dx, j 



3 > 

-1 



are known (note that, whenever m is even, this means that the first m — 1 Fourier coefficients of / are 
given. We will allow this minor discrepancy since it simplifies ensuing analysis). As a consequence 
of Theorem |2.4| we are free to choose the reconstruction space. The orthogonal projection of an 
analytic function onto the space P„_i of polynomials of degree less than n is known to converge 
exponentially fast at rate p~ n , where p > 1 is determined from the largest Bernstein ellipse within 
which / is analytic 9 . Hence, we let T„ = P n _i. Note that an orthonormal basis for T„ is given 
by the functions 

= \fI+\Pj{*), ieN, (3.1) 

where Pj is the j th Legendre polynomial. Moreover, if Q n is the orthogonal projection onto T„, 
then it is well-known that 

||/-Q„/|| <c f Vnp- n , (3.2) 

where c/ depends only on the maximal value of / on the Bernstein ellipse indexed by p. Naturally, 
we could also assume finite regularity of / throughout, with suitable adjustments made to the various 
error estimates. However, for simplicity we shall not do this. 



With this to hand, provided m > 8(n; 9) for some 9 > 0, where 0(n; 9) is defined in (2.22 1, the 



approximation f n ,m obtained from the reconstruction procedure satisfies ||/— / n , m || < c(9)\\f—Q n f\\ 



(see Theorem 2.4 ) . In particular, ||/ — f n .m\\ < c(9)cfy/np~ n . Hence, exponential convergence of 
fn.m- The key question remaining is how large m must be in comparison to n to ensure such rates of 
convergence. Resolving this question involves estimating the quantity C n>m , a task we next pursue. 

3.1 Estimates for 9 (n; 9) 

For both numerical and analytical estimates of Q(n;9) we need to select an appropriate basis of 
i-i (recall from Section [2] that Q(n;9) is independent of the basis used). A natural choice is the 



orthonormal basis (3.1 1 of scaled Legendre polynomials. Fortunately, in this case, the inner products 
{4>k,ij)j) (i.e. the entries of the matrix U) are known in closed form: 



(0 fc ,^) = (-i) fc y— M+§C7't), jez, keN, (3.3) 
where J m is the Bessel function of first kind. This follows directly from 

3m{z) = \(-i) m j e izx P m (x)dx, VzeC, (3.4) 
(see [TJ 10.1.14]), where j m is the spherical Bessel function of the first kind, given by 

3m(z) = y^J m+ i(z). 

With this to hand, we may compute C„ jln (and, in turn, 0(n; 9)) via the expression given in Lemma 
In Figure 1 we display the functions 0(n; |) and &(n; |) against n. Immediately, we witness 



2.15 

quadratic growth of 0(n;9) with n, a result we verify in this section. In doing so, we shall also 
derive an upper bound for Q(n;9) in terms of n and 9. This gives an explicit, analytic condition 
for quasi-optimal recovery. Whilst such a bound is completely robust (in that it holds for all n), we 
notice from Figure [l] that 0(n; 9), when scaled by nT 2 , quickly converges to an asymptotic limit. In 
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practice it is wasteful to use a larger value of m than necessary (or, conversely, for fixed m a overly 
pessimistic value of n). Hence, in the second part of this section, we will also derive an asymptotic 
bound for Q(n; 9). 

We commence as follows: 



Lemma 3.1. Suppose that T„ ~ P n _i, S m is the space spanned by the first m Fourier basis functions 
and m > max{2, - n}. Then C n ., n satisfies 



C > 1 



4(tt - 2)n 2 



TT^LfJ-l)- 

Proof. From the definition of C n%m and the fact that {ipj} is an orthonormal basis we have 
1-Cn,m = l- inf (V m <f>, </>) = sup (<j) - V m <j), 4>) = sup \\<j) - V m (j)\\ 2 , 



11 = 1 



11=1 



where V m is the Fourier projection operator. It now follows that 1 — C n , m < Sfc=o ll^fc — ^'m'/'fcll 2 , 
where 4>k is given by (3.1). By Parseval's theorem and the expression (3.3), we find that 

U k -vM\ 2 = £ k -^\J kH {j*)\ 2 . 

IjI>LxJ J 

Now, using a known result for Bessel functions |32j . it can be shown that 



2fc + 1 



provided jtt > k + ^. Hence, for m > ^n, 



^ 2(2fc + l) 



1 



(fc+|) 2 



(3.5) 



Now, it was shown in [32] that 



j>m 



3\JP- 



1 < i arcsin — -r- , whenever m > c + i . This gi 



gives 



^ - fm^fe < - arcsin 

7T 



2fe + l 



(2LfJ-l)7T 
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and 



n-l 



arcsm 



1 — C n m ^ ~ / , 

k=0 

We estimate this sum by the integral of arcsini. We have 

2 

l-C_ m <2(2LfJ-l) 



2k + 1 



(2L?J"1)t 



arcsint dt. 



Now, arcsinx < (arcsinl)a; for < x < 1. It follows that = J arcsin t dt < F(l)x 2 . Com- 

puting this integral exactly, we arrive at F(l) = ~ — 1. Upon substituting x — ^y^-i)* > * ms 
completes the proof. □ 

This lemma confirms that it is sufficient for to to scale quadratically with n for quasi-optimal 
recovery. Using this result, we find that 



Theorem 3.2. Suppose that T„ and S m are as in Lemma 3.1 Then, for n > 2, <d(n;9) satisfies 



@{n;9) < 2 



1 2(?r-2) 



2 tt 2 (1 — 6») 



Vn € N. 



Proof. Suppose that to > {2, — n}. Then, by Lemma 3.1 C„ im > 6 if 

1 4(7r-2)n 2 > 



TT^LfJ-l) 



Rearranging, we find that 
2 



> 



4(vr - 2)n 2 

1 + — — 

tt 2 (1-6») 



TO > 2 



1 2(?r-2) 



2 7T 2 (1 



and the theorem is proved, provided the right-hand side exceeds max{2, ^0}. Since n > 2, the 
right-hand side is certainly greater than 2. Moreover, 

4(7r-2)n 2 8(?r-2)n 2n 

1 H — — > — — > — , 

7T 2 (1 — 9) 7T 7T 

as required. □ 

Using a similar approach, we are also able to obtain an asymptotic bound for 0(n; 9) that is 
sharper than if were to use Theorem 3.2 directly: 

Theorem 3.3. Suppose that T n and S m are as in Lemma \3.1\ Then the function Q(n;9) satisfies 

4 



n- 2 Q(n;9) < 



7T 2 (1 - 



■O In- 



n—too. 



Proof. Suppose that to = cn 2 and recall (p.5L Since j < n and k > |cn 2 , we deduce that 



<p A ||=2 <r 2 ( 2fc + 1 ) V- 1 
»fc - "mMl < Z3 2^ ^2 



o 



_ 4 . _ 4(2fc + 1) 



7T- 



j>L^J 



+ G(n- 4 ). 



Hence 



A — /I 

1 - C n . c1l 2 < V(2fc + 1) + O (n- 2 ) = ^+0 (n- 2 ) 

fe=0 



Rearranging now gives the result. 



□ 
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In Figure [2] we compare the function n~ 2 0(n;6) for = 5, | and the global and asymptotic 
bounds of Theorems |3.2| and |3.3| Note that both bounds are reasonably sharp in comparison to 
the computed values. In particular, as n — > 00, n~ 2 0(n; |) quickly approaches the limiting value 
c s» 0.38, whereas the global and asymptotic upper bounds are 0.93 and 0.81 respectively. 

At this moment, we reiterate an important point. Whilst Legendre polynomials were used in 
the proof of Lemma |3.1[ the constant C„ jm is independent of the particular reconstruction basis, 
and is only determined by the space T„. Hence, Theorems |3.2| and |3.3| provide a priori estimates 
regardless of the particular implementation of the reconstruction procedure. In the next section, we 
discuss the choice of polynomial basis and its effect on the numerical method. 

Remark 5 In some applications, medical imaging, for example, oversampling is common. Formally 
speaking, this is the situation where we wish to recover a function / with support in [—1, 1] from its 
Fourier samples taken over an extended interval K D [—1, 1] (e.g. K = [—-,-] for some < e < 1). 
In this case, proceeding in a similar manner to before, we let H = L 2 (K), 



■tpj(x) 



x G K, 



where c = \\K\ and T„ = {</> : </>|[-i,i] € Pn-i, supp((/>) C [—1,1]}. Using similar arguments to 
those of Lemma 



and 



9(n;0) < 2 



3.1 one can also derive estimates for C n . m and 0(n; 9) in this case. In fact, 

4(tt - 2)n 2 

CK z (m — \) 

, VneN, n~ 2 Q(n;9) < 



(3.6) 



1 2(7T-2) 2 

2 ctt 2 {\-9) 



C7T 2 (1 — ' 



+ 0(n 2 ) , n -> 00. 



In particular, we retain the scaling m = O (n 2 ), regardless of the of size of the interval K. 



3.2 Choice of polynomial basis 

The results proved in this section are independent of the polynomial basis used for implementation. 
In selecting such a basis, there are two questions which must be resolved. First, how stable is the 



resultant method, and second, how can the entries of the matrix U (as defined in Section 2.3 ) be 
computed? A straightforward choice is the orthogonal basis of Legendre polynomials (3.1|. In this 
case, A — I, where A is the Gram matrix for {4>o, . . . , </>„_i}, making the method well-conditioned 
(Lemma 2.13). Moreover, the entries of U are known explicitly via (3.3|. 
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Having said this, there is also interest in reconstructing in other polynomial bases. In many 
circumstances it may be advantageous to have an approximation f n ,m that is easily manipulable. In 
this sense, an approximant composed of Legendre polynomials is not as convenient as one consisting 
of Chebyshev polynomials (of the first or second kind); the latter being easy to manipulate via the 
Fast Fourier Transform. To this end, the purpose of this section is to detail the implementation of 
this method in terms of general Gegenbauer polynomials. 

Gegenbauer polynomials arise as orthogonal polynomials with respect to the inner product 

(f,g) x = J /(z)sR(i - z 2 ) A -s dz, \>- 1 -. 

For given A, we denote the j th such polynomial by C A € Pj. Important special cases are the Legendre 
polynomials (A = |), and Chebyshev polynomials of the first (A = 0) and second (A = 1) kind. By 
convention [7] (see also [3T]), each polynomial C A is normalised so that 

A r(j + 2A) 

C i {l) - j\T(2X) ' (3 - ?) 
where T is the Gamma function, in which case it is known that (see [71 p. 174]) 



r(j + 2A)r(A + i) 

^•|lA-^ i!r(2A)r(A)(i + A) , (3.8) 



where ||/||a = y/{f, f) x - With this to hand > 

we now define 

1 



cH Cj\ j -0,1,2,..., (3.9) 



A 



and seek to reconstruct / in this basis. 

Our first task is to compute the entries of the matrix U. For this, we need to compute integrals 
of the form 

I k {z) = J C£(x)e izx dx, k = 0,1,2,..., 

where zei Fortunately, such integrals obey a simple recurrence relation: 
Lemma 3.4. For z =/= 0, the integrals Ik(z) satisfy 

/ 2=2C A 1 , h(z) = 2iC*(l) - 2 , 

z z A 

h+l{z) = 2 ^t±A Ik{z) + h _ l{z) i c ' z + { -^ lz [c A +1 (i) - ctm , k = i, 2, . . . . 

When z = 0, we have 

7 (0) - 2C A (1), 4(0) - \ + {k { ^ [C fe A +1 (l) - C#_i(l)] , k = 1,2,.... 
Proof. Recall the identity (see [7j p. 176]) 

^ (x) = 217Ta)^ [c ^ 1 ~ c '^ i] ' J' = 1 ' 2 '---- 

Substituting this into the expression for Iu{z) and integrating by parts gives 

W = [(C£ +1 (x) - CjU*)) ef]^ - [I k+1 (z) - I k ^{z)] ■ 
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Rearranging now gives the general recurrence for fc > 1. For k = 0,1, we merely note that Cq (x) 
C* A (1), C x {x) = C£(l)x and that 



1 i^r i „sinz f 1 irT sin z — z cos 2 
e dx = 2 , / xe dx = 2 s . 



The result for z = is derived in a similar manner. □ 

Using this recurrence formula, the matrix U can be formed in O (mn) operations. With this to 
hand, we now turn our attention to the condition number of A: 



Theorem 3.5. Let A be Gram matrix for the vectors {<fio, ■ ■ ■ , 4> n -i}> where <pj is given by (3.9). 
Then, n(A) — O (rJ 2A_1 l) as n — > oo. In particular, whenever </>o, . . . , 4>n-i arise from Chebyshev 
polynomials (of the first or second kinds), then k(A) — O (n). 

To prove this theorem, we first require the following two lemmas. For convenience, we will write 
L^(— 1, 1), A > — ~, for the space of square-integrable functions with respect to the Gegenbauer 
weight function (1 — £ 2 ) A ~5. 

Lemma 3.6. Suppose that — | < A < ~. Then, for all g 6 L°°(— 1, 1), we have \\g\\ < \\g\\\ and, for 
some c\ > independent of g, 

II.9IU < C ANi A+ ^Niir A - (3.io) 

Conversely, if A > i then, for all g E L°°(— 1, 1), \\g\\ < \\g\\\ and 

1 * 2. 

\\g\\<c x \\g\\f^\\g\\^. (3.11) 

Proof. Suppose first that — | < A < \. Trivially, ||<7|| < \\g\\\. Now consider the other inequality. 
For any < e < 1, we have 

NIa = J\(x)\ 2 {i-^) x -idx 

\g(x)\ 2 (l~x 2 ) x -i dx+ I \g{x)\ 2 {l-x 2 ) x -^dx 
<(l-(l- e ) 2 ) A -5|| g || 2 + 2|| 5 || 2 f (l-x 2 ) x 'idx, 

Jl-e 

where || • ||oo is the uniform norm on [-1,1]. Note that (1 - (1 - y) 2 ) A "5 < y A "5, Vy e (0,1). It 
now follows that 

\\9\\l<e x -i\\g\\ 2 + -^ T e x +i\\gf 0O , < e < 1. 

I II 2 

Let c > 2 be arbitrary. Then ||g|| 2 < cH^H;^, so we may let e = J^ 2 . Substituting this into the 
previous expression immediately gives (3.10|. 

Now suppose that A > |. Once more, trivial arguments give that \\g\\\ < ||<?||. For the other 
inequality, we proceed in a similar manner. We have ||g|| 2 < e2~ A || 5 || 2 + 2e||<^|| 2 . For c> 2 we now 

2 I I 

set e = ' which § ives (3 - 11 D 

Lemma 3.7. Let A > A. Then, for all <f> £ P„_i, Halloo < ^aII^Ha) where k n \ depends only on n 
and A and satisfies fc„,A = O ^ra A+3 ^ as n —> oo. 
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Proof. Let {<fio, . . . , (j> n -i\ be given by (3.9|, and write an arbitrary cf> € P n -i as </> = Y^ij=o a j j 
where Y^=o \ a j\ 2 ~ 11011a- t ne Cauchy-Schwarz inequality, 

n-l 

H\\l<H\\lY, UiWl = ^110111 
J=0 



We now wish to estimate ||0j||oo- Recall that ||C*||oo = Cf{\) 7, p.206]. Hence, by < 3.7 ) and (3.8 1, 
we have 

u ii 2 = r (j + 2A )(j + A ) 
ll<Pjllo ° j!^r(2A)r(A + 1)' 

Consider the ratio j ! 2A - ) . By Stirling's formula, 

Hence ll^j ll 2 ^ = (j 2A )> which gives fc 2 A = (n 2A+1 ), as required. □ 



Proof of Theorem \3.5\ Since A is Hermitian and positive definite, its condition number is the ratio 
of its maximum and minimum eigenvalues. By a simple argument, we find that 

Amax(^) = sup jnfr, A min (A) = inf f^lrr- 



Consider the case A > \. By Lemma 3.6 we have A m i n (j4) > 1 and 



A m a X (^4) < c 2 sup 



2^4 



Using Lemma |3.7| we deduce that A max (A) = (n 2 * -1 ), as required. For the case —\ < A < |, we 
proceed in a similar manner. □ 

This theorem confirms that the method can be implemented using Chebyshev polynomials whilst 
incurring only a mild growth in the condition number. It follows that, if conjugate gradients are 
used to compute the approximation, the total computational cost of forming f n ^ m is 0(mn5), as 
opposed to O (mn) in the Legendre polynomial case. In the next section we present several examples 
of this implementation. 



Remark 6 Whilst Theorem 



3.5 



provides an asymptotic estimate for n(A) (and hence k(A)), it 
may also be useful to derive global bounds. With effort, one could obtain versions of Lemmas |3.6| 
and |3 . 7| involving explicit bounds. For the sake of brevity, we shall not do this. However, whenever 
Chebyshev polynomials are used (arguably the most important case), it is a relatively simple exercise 
to confirm that 

k(A) < 2\/2n, k(A) < 3371--3 [n(n + |)(n + 1)] 3 , 
in the first and second kind cases respectively. 
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Figure 3: Error in approximating f(x) — e x cos 4a; by f n ,m(x) for n = 1, . . . , 40. Left: log error log 10 ||/ — 
/n,m||oo (squares) and log 10 ||/ — /»,m|| (circles) against n. Right: log error against m = 0.2n 2 . 




Figure 4: Error in approximating f(x) = e x cos 4x by f n ,m(x) (squares) and Q n f{x) (circles) for n = 
1, . . . , 40. Left: log uniform error. Right log L 2 error. 



3.3 Numerical examples 

We now present several numerical examples of this method. All examples employ the value m = 
0.2n 2 , and the first series of examples consider the implementation using Legendre polynomials. In 
Figure [3] we consider the function f(x) = e~ x cos4a;. Since / is analytic in this case, we witness 
exponential convergence in terms of n and root exponential convergence in terms of m. Note the 
effectiveness of the method: using less than 100 Fourier coefficients, we obtain an approximation 
with 13 digits of accuracy. 

As indicated by Theore m |2.4[ the approximation f n . m is quasi-optimal. To highlight this fea- 
ture of the method, Figure |4| displays both the error in approximating / by f n , m and the best 
approximation Q n f- Note the very close correspondence of the two graphs. 

The example in Figures [3] and [4] is, in fact, entire. Hence, the approximation f n . m converges 
super-geometrically in n (as seen in Figure[3]). For a meromorphic function, with complex singularity 
lying outside [—1, 1], the convergence rate is truly exponential at a rate p. This is demonstrated in 
Figure jHJ the approximated function being f(x) = j^s- Note that, despite the poles at x = ±i, the 
approximation f n . m still obtains 13 digits of accuracy using only 250 Fourier coefficients. 

Next we consider reconstructions in other polynomials using the work of the previous section. 
In Table [l] we give the error in approximating the function f(x) = e x cos Ax with Chebyshev poly- 
nomials of the first and second kinds. Note that the resulting uniform error is virtually identical to 
the case of the Legendre polynomial implementation (unsurprisingly, since all three implementations 
compute exactly the same approximation f n , m , up to numerical error). Moreover, as evidenced by 
Table pi the payoff is only mild growth in the condition number k(A). 
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n 


5 


10 


15 


20 


25 


30 


35 


40 


(a) 
(b) 
(c) 


1.45e0 
1.45e0 
1.45e0 


1.85e-3 
1.85e-3 
1.85e-3 


3.03e-7 
3.03e-7 
3.03e-7 


2.53e-12 
2.53e-12 
2.49e-12 


1.06e-14 
3.51e-14 
6.76e-14 


8.42e-14 
1.16e-13 
7.33e-14 


4.06e-14 
4.57e-14 
6.40e-14 


5.31e-14 
7.70e-14 
5.15e-14 



Table 1: Comparison of the error ||/ — f n ,m\\oo with m = 0.2n 2 , where f„ y m is formed from (a) Legendre 
polynomials and Chebyshev polynomials of the (b) first and (c) second kinds. 



3.4 Connections to earlier work 

Rather than choosing m such that C n ^ m > 6, it may appear advantageous to find the minimum m 
such that C„. m > 0. In other words, the smallest m such that /„ jm is guaranteed to exist. Letting 



9 = in Theorems 3.2 and 3.3 we immediately obtain a sufficient condition of the form m > en 2 , 
for some c > 0. However, this result is far too pessimistic: it is known that reconstruction is always 
possible, provided m > n [32; . For this reason, it may appear favourable to reconstruct using 
to = n, leading to the so-called inverse polynomial reconstruction method |37[ 138], Unfortunately, 
however, this approach is extremely unstable. The linear system has geometrically large condition 
number, making the procedure extremely sensitive to both noise and round-off error. Moreover, a 
continuous analogue of the Runge phenomenon occurs. Roughly speaking, the approximation f n . m 
only converges to / if geometric decay of ||/ — Q n f\\ is faster than the geometric growth of ||.A _1 ||, 
meaning that only functions analytic in sufficiently large complex regions can be approximated 
by this procedure (as discussed in detail in [3j, this behaviour can be understood in terms of the 
operator-theoretic properties of finite sections of certain non-Hermitian infinite matrices). On the 
other hand, by allowing m to range independently of 77, we overcome all these difficulties, and obtain 
a stable method whose convergence is completely determined by the convergence of Q n f to /. 

For the specific example of Legendre polynomial reconstructions from Fourier samples, this ap- 
proach has also been recently considered in [32] • Therein, the estimate m = O (n 2 ) was derived, 
along with bounds for the error. Naturally, this problem is just one specific example of our general 
framework. However, within this context, our work improves and extends the results of 32 in the 
following ways: 

1. Reconstruction is completely independent of the particular polynomial basis used. In partic- 
ular, the estimates for G(n; 0) and ||/ — / n , m || are determined only by the spaces T n and S m . 
This allows for analysis of reconstructions in arbitrary polynomial bases, not just the Legendre 
polynomials used in |32j . 



The estimates for 0(n;#) in Theorems 3.2 and 3.3 improve those given in |32| . In particular, 
it was shown in Theorem 4.2] that 

8 1 

C nan 2 > 1 arcsin — , Vn € N, a > 1, (3.12) 

7r net 
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n 


5 


10 


15 


20 


25 


30 


35 


40 


(a) 


3.57 


5.55 


4.21 


5.20 


4.40 


5.06 


4.50 


6.77 


(b) 


13.74 


49.99 


52.63 


91.89 


92.89 


133.02 


133.49 


191.19 


(c) 


3.90 


5.67 


7.25 


9.33 


11.91 


13.96 


16.56 


18.92 



Table 2: Comparison of the condition number k(A) with m = 0.2n 2 , where A is formed from (a) Legendre 
polynomials and Chebyshev polynomials of the (b) first and (c) second kinds. 



(our constant C„, m corresponds to the quantity c 2 m in [32]). Conversely, Lemma 3.1 leads to 
the improved bounds 



4(tt - 2) 
7r 2 (o! — nr 



Vn > max 




a > 0, 



and 



,2 > 1 



01 



OO, 



a > 0. 



(3.13) 



(3.14) 



Not only are these bounds sharper, they also hold for a greater range of a, thus permitting 



reconstruction with m 



for any a > 0, as opposed to just a > 1. This leads to savings 



in computational cost, and, in cases where m is fixed, allows larger values of n to be used, 



thereby increasing accuracy. To illustrate this improvement, note that (3.12 ) gives the estimate 
k(A) < 5.71 when m ~ n 2 . Conversely, our estimate (3.13) yields the bound 2.61 for n > 2, 



and (3.14) gives the asymptotic bound 1.68. To compare, direct computation of k(A) indicates 
that k(A) < 1.32 for all n, and k(A) — > 1.2 as m — > oo. 

Piecewise analytic functions and function of arbitrary numbers of variables can be recovered 
in a analogous fashion, with similar analysis (see Sections 3.5 and [4] respectively). 



3.5 Reconstruction of piecewise analytic functions 

Naturally, whenever the approximated function is not analytic, the convergence rate of the polyno- 
mial approximant f n , m to / is not exponential. For example, consider the function 

f(x) _ / (2eM*+D ~ 1 - e*)(e* - l)" 1 x e [-1, -A) , 

This function was put forth in |46j to test algorithms for overcoming the Gibbs phenomenon. Aside 
from the discontinuity, its sharp peak makes it a challenging function to reconstruct accurately. 
Since this function is discontinuous, we expect only low-order, algebraic convergence of f n ^ m in the 
L 2 norm, but no uniform convergence, an observation confirmed in Figure [6j 

However, by reconstructing this function in a polynomial basis, we are not exploiting the known 
information about /: namely, the jump discontinuity at x = — h. The general procedure set out in 
Section [2] allows us to use such information in designing a reconstruction basis. Naturally, since / is 
analytic in the subintervals [— 1, — |] and [— a better choice is to reconstruct / in a piecewise 
polynomial basis. The aim of this section is to describe this procedure. 

Seeking generality, suppose that / : [— 1, 1] — > R is piecewise analytic with jump discontinuities 
at — 1 < X\ < . . . < xi < 1. Let xq = — 1 and xi+i = 1. We assume that / has been sampled via 
fj = (/, ipj), j = 1, . . . , m, where (•, •) is the Euclidean inner product on L 2 (— 1, 1). In examples, 
these will be the Fourier samples of /, but the construction described below holds for arbitrary 
sampling bases consisting of functions defined on [—1, 1]. 

Throughout we shall assume that the discontinuity locations x\, . . . ,xi are known exactly. Al- 
though this may be a reasonable assumption in many applications, a fully-automated algorithm 
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must also incorporate a scheme for singularity detection. We shall not discuss possible approaches 
for doing this, and we refer the reader to 45j for further details. 

Given the additional information about the location of the singularities of /, we now design a 
reconstruction basis to mirror this feature. We shall construct such a basis via local co-ordinate 
mappings. To this end, let I r — [a; r ,x r +i], c r = h(x r +i — x r ) and define A r (x) — x ~ Xr ~ 1, so that 
A(/ r ) = [—1, 1]. Suppose now that T' n is a space of functions defined on [—1, 1] (e.g. the polynomial 
space P„_i). By convention, we assume that each <j> £ T' n is extended by zero to the whole real 
line, i.e. 4>(x) = for x £ 1, 1]. Let T„, r be the space of functions defined on I r , given by 
T„ r = {(f> o A r : <p £ T^}. We now define the new reconstruction space in the obvious manner: 

i 

T„ = {<j> : <f>\ Ir £ T Ur ,r, r = 0,...,l}, n = ^ n r , 

r=0 

and seek an approximation f n . m £ T„ to / via the conditions a m (f n>m ,4>) = a m (f,<j>), V0 £ T 



n ' 



where a m is defined in (2.5). Suppose now that {(pi, ...,</>„} is a collection of linearly independent 
reconstruction functions with T' n — spanj^i, . . . ,4> n }. We construct a basis for T„ by scaling. To 
this end, we let cf> r j = -^=4>j ° A r , and notice that T„ = span {(f> r .j : j = 1, . . . , n r , r = 0, . . . , I}. 
Note that, if {<fij} are orthonormal, then so are {<f> r ,j}- With this basis in hand, the approximation 

fn,m 

is now given by 

l n r 

fn.m — ^ ^ ^rjtftrj ■> 
r=0 j = l 

where the coefficients a r .j are determined by the aforementioned equations. As before, this is 
equivalent to the least squares problem Ua « / with block matrix U = [t/i, . . . , U{\, where U r is the 
m x n r matrix with (j, k) th entry 



1 



x r+1 



3 \ 



Here / = (/i, . . . , / m ) T , a = [a , ...,a t ] and a r = (a r ,i, . . . ,a r ,„J T . 

Naturally, estimation of the constant C n _ m is vital. The following lemma aids in this task: 

Lemma 3.8. The constant C n _ m satisfies 

i 

r=0 

where C„ r>TOir = inf^ eT „ r (P m <f>, 4>) ■ 

ll*ll=i 
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Proof. For g T n , denote 0|/ r by 0M. Assume that is extended to [—1,1] by zero, so that 
4> = Et= Since <^ W - 1 [s] for r ^ s > it; follows that 



1 - Cn,m = sup • 



EU{^-W [r W [rl >,, [r] 



£T W , r = 0,...,l, £||4>H||V0 



r=0 



Note that, for a r > and 6 r > 0, r = 0, . . . , I, the inequality 

i i i 



r=0 r=0 



holds. Setting a r = (^M — V m <t^ r \ 4>^) anc l = ||0' r '|| 2 and using this inequality gives 



i ^ / /v-<^ [rl -^ W ,^ r] ) 
1 - C n . m < sup < 2^ 



||0W||2 



G T- 



o,...,z, 



sW||2 



r=0 



< ^ sup j - 



^T nr , r , \\<t>\\ ^0 



and this is precisely ^ r=0 (l — C nr<n 



□ 



Let us now focus on piecewise polynomial reconstructions from Fourier samples, in which case T„ 
is the space of functions 4> with ip\j r a polynomial of degree n r . Regarding the rate of convergence 
of the resulting approximation f nm , it is a simple exercise to confirm that 



\\f-fn,m\\<c(9)CfJ2 



n r p r 



where c(9) is defined in (2.23), c/ is a constant depending on / only and p r is determined by the 
largest Bernstein ellipse (appropriately scaled) within which the function f\j r is analytic. Hence, 
exponential convergence of f n ,m to /. The main question remaining is that of estimating the function 
0(n; 9) for this reconstruction procedure. For this, we have the following result, which extends 
Theorems 13.21 and 13.31 to this case: 



Theorem 3.9. The function 0(n; 9) satisfies 



6(ra;0) < 2 



1 2(7T -2) 



E? 



2 n 2 (l-9)^cr 



Vn = 2jn r , no,..., ni e N, 



and 



Q(n;9) < 



1 2 



7r 2 (l 



E 



■0(1) 



r=0 



"0) 



, n; — ?> oo. 



Proof. In view of Lemma 3.8 it suffices to consider C nrtmrtr . To this end, let J = [a,fj] C [— 1, 1], 
T n ,j be the space of functions <p with supp(0) C J and 0| j € P n -i, and define 



n,m Sup 



Let A(x) = — 1, where c — |(/3 — a), and write ^ = $oA, where supp(<£>) C [—1, 1]. Consider 
the quantity (4>,ipj). By definition of V'j, we have 



{<i>Ai) 



'2 J- 



0(x)e- ij7ra: Ax 



7=2 



-ij7r(a+c) 



A(l) 



$(y)e" 



i]Trcy 



dy. 



A(-i) 



25 




500 





-B.5 




0.5 


L0 


























-16 







Figure 7: Error in approximating the function (3.15 1 by / n , m (x). Left: log error log 10 ||/— /„,. 
and 



(squares) 

/ 2 2 \ 

Jio 11/ ~ fn.m || (circles) against m = 1, . . . , 500 with no = ni chosen so that m = | ( ^ + ) • Right: 
the error log 10 |/(x) — f n , m (x)\ against a; € [— 1, 1] for m = 20, 40, 80, 160. 



Let K = [A(— 1), A(l)] = A([— 1, 1]) 2 [—1, 1] and let V m ,K be the Fourier projection operator based 
on the interval K. It now follows that 

E„ >m = sup{<$--p ro ,K$,$) :supp($) C [-1,1], $![_!,!] €P„_i, ||*|| = 1}. 

This is now precisely the setup of Remark [5] Using (3.6), we therefore deduce that 

4(tt - 2)n 2 
"' mS C7T 2 (2LfJ - I)' 

Letting J — I r , c — c r and using Lemma |3.8| we now obtain 



C > 1 



' 2 



4(tt - 2 
7r 2 (2Lf J - 1) < 



from which the result follows immediately. 



(3.16) 
□ 



To implement this scheme, it is necessary to compute the values (3.5). By changing variables, it 
is easily seen that 



where d r — h(x r +i + x r ). Since (3.4) holds for all z G C, it follows that 



-ii7Td r (_ i )fcJ fc + 2 



(3.17) 



whenever the functions <j) r ^ arise from scaled Legendre polynomials. Naturally, if the functions 



arise from arbitrary scaled Gegenbauer polynomials, computation of the values (3.5) can be carried 



>r,k 



out recursively via the algorithm described in Section 3.2 (for appropriate choices of z). 

In Figure [7] we apply this method to the function (3.15) using the orthonormal basis of scaled 
Legendre polynomials. The improvement over Figure |6| is dramatic: using only m w 250 (with 
uq = rix = 16) we obtain 13 digits of accuracy. Note that, as expected, root exponential convergence 
occurs. Moreover, as predicted by (3.16), and illustrated in Table [3j the condition number of the 
matrix A remains small. 
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m 


10 


20 


40 


80 


160 


320 


640 


1280 


k(A) 


0.05 
19.88 


0.34 
2.92 


0.33 
3.06 


0.44 
2.27 


0.44 
2.27 


0.47 
2.11 


0.49 
2.03 


0.50 
1.98 



Table 3: The constant C n ,m an d the condition number k(A) against m. 



3.6 Comparison to existing methods 

Numerous methods exist for recovering functions to high accuracy from their Fourier data. Appli- 
cations are myriad, and range from medical imaging [5J [5] to postprocessing of numerical solutions 
of hyperbolic PDEs [UJ ES] • Prominent examples which deliver high global accuracy (in contrast to 
standard filtering techniques, which only yield high accuracy away from the singularities of a func- 
tion [3D]) include Gegenbauer reconstruction [551 EM ED , techniques based on extrapolation [T9"] . 
Pade methods [18] and Fourier extension methods [101 [33] (for a more complete list, see [12] and 
references therein). 

Whilst many of these methods deliver exponential convergence in terms of to (the number of 
Fourier coefficients), they all suffer from ill-conditioning. This comes as no surprise: the problem of 
reconstructing a function from its Fourier coefficients can be viewed as a continuous analogue of the 
recovery of a function from m equidistant pointwise values. As proved in [42], any method for this 
problem that converges exponentially fast in m will suffer from exponentially poor conditioning. We 
conjecture that a similar result holds in the continuous case. 

Aside from increased susceptibility to round-off error and noise, ill-conditioning often makes 
so-called inverse methods (e.g. extrapolation and Fourier extension methods) costly to implement. 
Conversely, the method proposed in the paper does not suffer from any ill-conditioning. This negative 
consequence of |42) is circumvented, precisely because we witness only root exponential convergence 
in m. However, an advantage of this approach is that it delivers exponential convergence in n, 
the degree of the final approximant f n , m - In many applications it may be necessary manipulate 
fn,mi its relatively low degree making such operations reasonably cheap. Thus, this method has the 
advantage of compression, a feature not shared by the majority of the other methods mentioned 
previously. 

A well-established and widely used alternative to this method is the Gegenbauer reconstruction 
technique, proposed by Gottlieb et al [3T]. Much like this approach, it computes a polynomial 
approximation. Yet it stands out as being direct, meaning that no linear system or least squares 
problem is required to be solved. Whilst the original method has been shown to suffer from a gen- 
eralised Runge phenomenon .11], thereby severely affecting its applicability, an improved approach 
based on Freud polynomials was recently developed in [24] , 

Comparatively speaking, this approach delivers exponential convergence in O (m 2 ) operations. 

On the other hand, our method obtains root exponential convergence at a cost of 0(rrfi) oper- 
ations. However, despite being theoretically more efficient, the various constants involved in the 
Gegenbauer/Freud procedure tend to be rather large. 

Indeed, in Table g we compare the error in approximating the function (|3.15[) using both pro- 



cedures (the data for the Freud method is taken from [Ml Table 1]. Note that the parameter N 
used therein is such that 2N = to is the total number of Fourier coefficients). As is evident, the 
method proposed in this paper obtains an error of order 10 -14 using less than 256 Fourier coeffi- 
cients, whereas the Freud procedure does not reach this value until more than 1024 coefficients are 
used. 

The most likely reason for this improvement is that the generalised reconstruction method is 
quasi-optimal, thereby delivering a near-optimal polynomial approximation, whereas the Gegenbauer 
and Freud procedures do not possess this feature. Indeed, in the Gegenbauer case at least, we can 
quantitatively explain this phenomenon (the corresponding results for the Freud polynomial case 
are not so well understood). It is known that, for analytic functions /, the approximation F„ jm 
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m 


64 


128 


256 


512 


1024 


2048 


4096 


(a) 
(b) 


8.90c-01 
2.40e-04 


1.37e-01 
8.36e-09 


1.84e-04 
2.40e-14 


1.01e-07 
1.38e-14 


9.33e-13 
1.74e-14 


5.27e-13 
2. 26c- 14 


5.23e-14 
2.59e-14 



Table 4: Comparison of the (a) Freud polynomial method and (b) generalised reconstruction method applied 
to (3.151. Here m is the total number of Fourier samples. 



obtained from the Gcgcnbauer procedure converges geometrically at rate 7 1 ™, where 

_ t {fi + 2af+ 2a \ 

(see [301 eq n - (4-12)]). Here A = am is the parameter of the Gegenbauer polynomials used and 
n — /3m is the degree of the polynomial F n>m . Thus, in practice, the Gegenbauer method may 
converge much more slowly than the best polynomial approximation Q m f, which converges at rate 
p~ m . Conversely, our approach offers root exponential convergence at a rate determined solely by 
p: ||/-/„, m || ~p-^. 

Suppose now that the number of Fourier coefficients is fixed. Then, ignoring various constants, 
the Gegenbauer method will offer a lower error than the generalised reconstruction procedure only 
when 7™ < p _v ™. In other words, 




In typical examples (see [501 I2I])> the parameter values a = /3 = \ have been used, giving the 
condition m > 19(logp) 2 . Thus, for reasonably analytic functions / (i.e. analytic in a sufficiently 
large Bernstein ellipse), the Gegenbauer procedure will only begin to outperform the generalised 
reconstruction method when the parameter m is rather large. Moreover, whenever / is entire (as 
is the case with the example in Table [4]), the generalised reconstruction procedure will converge 
super-geometrically (in n = 0{^/m)), whereas the Gegenbauer method may still converge only 
exponentially at rate 7. Thus, in this situation the Gegenbauer method may never outperform the 
generalised reconstruction method for any finite m. 

Aside from improved numerical performance, let us mention several other benefits. First, as 
discussed, the final approximation consists of only O(y'rn) terms, as opposed to O (m). Second, the 
basis for the polynomial reconstruction space T„ can be chosen arbitrarily (in particular, indepen- 
dently of m) without affecting the convergence. The only downside is a mild increase in condition 
number if nonorthogonal polynomials are employed. In contrast, for the Freud/ Gegenbauer proce- 
dure, only very specific types of polynomials can be used (which may not be simple to construct 
or manipulate 24j), and, whenever the number of samples m is varied, all polynomials used for 
reconstruction must be changed. 

To somewhat temper our claims, we mention that we have only considered one particular test 
function. There may well be problems where the Freud/Gegenbauer procedure performs better, and 
the intention of future work is to present a more thorough comparison of the two methods. That 
aside, one advantage of Gegenbauer method is that it is local: the approximation in each subdomain 
is computed separately and independently of any other subdomain. Conversely, with our approach, 
the computations are inherently coupled. Nevertheless, it may be possible to devise a local version 
of our approach, a question we intend to explore in future investigations. 



max 



4 Reconstructions in tensor-product spaces 

Thus far, we have focused on the reconstruction of univariate functions from Fourier samples. A 
simple extension of this approach, via tensor products, is to functions defined in cubes. The aim of 
this section is to detail this generalisation. 



28 



To formalise this idea, let us return to the general perspective of Section [2] Suppose that 
the Hilbert space H can be expressed as a tensor-product H = Hi ® • • • ® ELj of Hilbert spaces 
Hi, i — l,,..,d, each having inner product (•,•),. Note that, for / = /i ® ••• ® /d € H and 
.9 = <?i ® • • • <8> 5d G H, we have 

d 
i=l 

Now suppose that the sampling basis consists of tensor-product functions. To this end, let 

ipj = i>i, h <g> •• • ® Vd,i d , j = (ii> ■ ■■,jd)£ N d , 

and, for m = (mi, . . . , m<j) G N d , set 

«S m = span : j = (ji, • . . , jd), l<ji<m h i = l,...,d}. 

We assume throughout that the collection {V'jjl^Li is a Riesz basis for Hi for i = 1, . . . , d (hence 
{ipj } is a Riesz basis for H) . With this to hand, we define the operator V m : H — > S m by 

mi rcid 

ji=i jd=i 

Note that P m (/i ® • • • <8> /d) = ?i, mi /i <& • • ■ ® Vd,m d fd, where 7> ijm< : H^ -> S ijm< is defined in the 
obvious manner. In a similar fashion, we introduce the reconstruction vectors <j),j — 4>\,j 1 ®- ■ ■®4>d,j d i 
which form a basis for the reconstruction space 

T n = span{0j : j = (ji, . . .Jd), l< j l <n i , i = 1, . . . ,d} , n = (n x , . . . ,nd) £ N d . 

As before, we construct the approximation f n>m £ T„ via the equations a m (/„. m , </>) = a m (f,(j)), 
V0 e T n , where a m (/, g) = ("P m /, g), V/, g € H. 

To cast this problem in a form suitable for computations, let J/W 6 C"»iX"» k e the matrix with 
(j, fc) th entry (ipij, 4>i,k) ■• Let U £ C m '™ be the matrix of the d-variate reconstruction method, where 
to = toi . . . rrid and n = n\ . . . rid- Then, it is easily shown that 

d d 
i=l i=l 

where A = WU, and A® = (£/M)t[/M ; and, in this case, Bi <g> B 2 denotes the Kronecker product of 
matrices the B\ and B 2 . By a trivial argument, we conclude that the number of operations required 
to compute f n .m is of order (raimi) . . . (ridTUd) y/ k(A) . 

Recall that the spectrum of the Kronecker product matrix B\ ® B 2 consists precisely of the pairs 
A/i, where A is an eigenvalue of B\ and fi is an eigenvalue of B 2 . From this, we easily deduce that 

d 

k(A) = JJk(AW). 

i=l 

Hence, n(A) is completely determined by the matrices A^ , with the i th such matrix corresponding 
to the univariate reconstruction problem with sampling basis an d reconstruction basis 

{<^ij}"=i- Unsurprisingly, a similar observation also holds for the quantity C n ^ m : 

Lemma 4.1. The constant C n . m satisfies C n _ m = Jlf=i Ci,ni,mr 



Proof. By Lemma |2.15| C n . m = A^A^A) and C i)rH , mi = : A min ((iW)- 1 AH), 1 = 1,...^, where 
A and AM are defined in the obvious manner. Since A = AM ® ■ ■ ■ ® A^, the matrix A -1 is the 
Kronecker product of matrices (AM) -1 . The result now follows immediately. □ 
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4.1 Reconstruction of piecewise smooth functions 

Having presented the general case, we now turn our attention to the reconstruction of a piecewise 
smooth function / : [— 1,1] — > K. We shall make the significant (see later discussion) assumption 
that / is smooth in hyper-rectangular subregions of [—1, l] d . To this end, let li G N for i = 1, . . . , d 
and suppose that 

-1 = xo,% < xi t i < . . . < xi ui < x u+ i, :l = 1, 

and define l r:l = [x r .i, av+i.i] for r = 0, . . . , l{. For r = (ri, . . . , r d ), we write I r = I Tl .i X • • • x I rdld , 
so that the collection 

{I r : r = (n,...,rd), n=Q,...,k, i = l, ...,d], 

consists of disjoint sets whose union is [— 1, l] d . We assume that / is smooth within each subdomain 
I r . In addition, for r*j = 0, . . . , h and i = 1, . . . , d, let c ru i — 5(2:^+1,$ — x rtl i) and set A ri .i(x) — 
x ~ Xri,i - 1, x G I n . Note that A n ^(I rui ) = [-1, 1]. 

We now construct a reconstruction space. To this end, for n G N let T' i5 dimT^ = n, be a space 
of functions <f> : K —> C with supp(0) C [—1,1]. Define 

T w = {0oA r , t : cfreT'J, neN. 

Now suppose that n is the vector (n\, . . . , n^), where 

h 

71% — ^ Tl>r,i j ^ — 1 ) • ■ • j d) 

r=0 

for some n r ^ G N. We now define the reconstruction space T„ by 

d h 

We seek a basis for this space. Let {(pi, . . . , <fi n }, n € N, be a basis for T^, and set 

</>rj,i = _ <ftj A ri j. 

A basis for T„ is now given by 

{0ri ,j u l ® • • • ® <j>r d ,j d ,d, j = 1, • • • j = 0, . . . i = 1, . . . , d} . 

This framework gives a general means in which to construct reconstruction bases in the tensor- 
product case for functions which are piecewise smooth with discontinuities parallel to the co-ordinate 
axes. Suppose now that we consider the recovery of such a function from its Fourier samples. Using 
the above framework, we are easily able to construct a basis consisting of piecewise polynomials of 
several variables. The main question remaining is that of estimating the function C n>m . However, 
in view of the Lemma |4.1| and the results derived in Section |3.1| a simple argument gives 

Theorem 4.2. Suppose that n = (ni, . . . , rid), where rij = X)r=o n "> an( ^ ^ 

6(n;0) = min{m = (mi, . . . , m d ) : C„, m > 9}, 9 G (0,1). 

If 9\, . . . ,9d G (0, 1) satisfy 9 — 9\ . . . 9 d , then we may write 

G(n;9) = (0i(ni; 9 X ), . . . , @ d (n d ; 9 d )) , 
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where, for i = 1, . . . , d, 



and 



Q l {n l -9 l ) < 2 



1 2(tt - 2) 



2 71-2(1-6),;) 



E 



4 n 2 . , 

©i(n»; 0j) < 2 , , 2^ — - + O (1) , n ,i, . . . -> oo. 

^ V 1 "ij „_ n ^i* 



The main consequence of this theorem is the following: regardless of the dimension, the variables 
mi,..., must scale quadratically with ni,...,n<j to ensure quasi-optimal recovery. Consider 
now the most simple example of this approach: namely, where / is smooth in [— 1, l] d , so that T„ 
consists of multivariate polynomials. In Figure [8] we plot the error in approximating the functions 
f(x,y) = e x y and f{x,y) — smSxy, using parameters mj — m 2 = 0.5nf and = n±. As in the 
univariate case, we observe the accuracy of this technique. For example, using only mi = 7712 ~ 200 
and Tii = 7i2 ~ 20 we obtain an error of order 10~ 14 . 

Remark 7 This approach (and many others based on tensor-product formulations) has the signif- 
icant shortcoming that it requires the function to be singular in regions parallel to the co-ordinate 
axes. Naturally, this is a rather restrictive condition. For a function with singularities lying on a 
curve (in two dimensions, for example), one potential alternative is to apply the one-dimensional 
method along horizontal and vertical slices, and recover the two-dimensional function from the 
resulting one-dimensional reconstructions. 

However, the generality of the reconstruction framework presented in this paper allows for another 
approach. Given that we can reconstruct in any basis, an obvious alternative is to subdivide the 
domain into triangular elements and use a finite element basis for reconstruction. This is a subject 
for future investigation. 



5 Other sampling problems 

Overcoming the Gibbs phenomenon in Fourier series is an obvious application of the general frame- 
work developed in Section [2j However, there is no reason to restrict to this case, and this framework 
can be readily applied to design effective methods for a variety of other problems. In this section 
we describe several related problems, and the application of this framework therein. 
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5.1 Modified Fourier sampling 

Modified Fourier series were proposed in [34j as a minor adjustment of Fourier series. In the domain 
[— 1, 1], rather than expanding a function / in the classical Fourier basis 

{cos jVa; : j £ N} U {sinjTvx : j G N + }, 

we construct the modified Fourier expansion using the basis 

{cos jnx : j e N} U {sin(j — ^)ttx : j e N + }, 

instead. Though this basis arises from only a minor adjustment of the Fourier basis, the result is an 
improved approximation: the modified Fourier series of a smooth, nonperiodic function converges 
uniformly at a rate of O (m _1 ), whilst Fourier series suffers from the Gibbs phenomenon. Although 
the convergence rate remains slow, the improvement over Fourier series, whilst retaining many of 
their principal benefits, has given rise to a number of applications of such expansions. For a more 
detailed survey, we refer the reader to [4]. 

We shall consider modified Fourier expansions in a somewhat different context. Given the sim- 
ilarity between the two bases, it is reasonable to assume that any sampling procedure (e.g. an 
MRI scanner) can be adapted to compute the modified Fourier coefficients of a given function (or 
image/signal), as opposed to the standard Fourier samples. Indeed, if 

JF/(f) = J /(aOe-^dc, 

is the Fourier transform of /, then the modified Fourier coefficients are precisely 

f 1 1 

ff = J f(x) cos jTnr dx = - [Ff(j) + F f{-])\ , 

// = J f(x) sin(j - ±)irxdx = - [Tf(j - §) + Ff{\ - j)] , 

and hence can be computed from samples of the Fourier transform. This raises the question: given 
that the general framework can handle sampling in either, is there an advantage gained from sampling 
in the modified Fourier basis, as opposed to the Fourier basis? As we shall show, provided the 
function is analytic and nonperiodic, this is indeed the case. Specifically, when we reconstruct in a 
polynomial basis, we require fewer samples to obtain quasi-optimal recovery to within a prescribed 
degree. 

Suppose that we carry out the reconstruction procedure as in Section|3]but using modified Fourier 
samples instead of Fourier samples. For this, we set 

1 A Lf J 

V m f{x) = -/ C + [ff cos jnx + ff sinjTrx . 

3=1 

Naturally, we consider the function 0(n; 9) once more. In Figure|9]we plot the function 0(n; 9) for the 
modified Fourier basis. Upon comparison with Figure [2j we conclude the following: using modified 
Fourier sampling, as opposed to standard Fourier sampling leads to a noticeable improvement. 
Specifically, n~ 2 Q(n; |) is approximately 0.15 for large n in the modified Fourier case, as opposed 
to 0.4 in the Fourier case. 

This result means that, if the number of samples m is fixed, we are able to take a much larger 
value of n in the modified Fourier case, whilst retaining quasi-optimal recovery (with constant c{9)). 
As a result, we obtain better, guaranteed accuracy. To illustrate this improvement, in Figure [10] 
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we compare the errors in approximating the function f(x) = e _a: cos8a: from either its Fourier or 
modified Fourier data. In both cases the number of samples m was fixed, and n was chosen so that 
the parameter C n ,m > \. As is evident, the method based on modified Fourier samples greatly 
outperforms the other. For example, using only m — 120 samples, we obtain an error of order 10 -14 
for the former, in comparison to only 1CP 4 for the latter. 

As in the Fourier case, to implement the modified Fourier-based approach it is necessary to have 
estimates for the function Q(n;9). These are particularly simple to derive: 

Lemma 5.1. Let T„ = P„_i and S m be the space spanned by the first m modified Fourier basis 
functions. Then the function 0(n; 6) satisfies 



9(n; 



< 



where k„ = 



sup 



Proof. In [5] it was shown that \\4> — Vm^W < 
result now follows immediately from the definition of C n 



2 

m7r 1 



for all sufficiently smooth functions 



The 

□ 



As a result of this lemma, analytical bounds for 6(n; 9) are dependent solely on the constant k n 
of the Markov inequality \\(f>'\\ < /c„n 2 ||(/)||, \/cf) E P„_i. Markov inequalities and their constants are 
well understood. The question of determining k n was first studied by Schmidt 43 , in which the 
estimates 

k n < — p, Vn, K n — >■ — , n — > co, (5-1) 

V2 7T 
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were derived. In |44j the following improved asymptotic estimate was also obtained: 



(n - 



1^2 r 

2 



1 - 



- 3 



12(n+i)a 



n > 5, 



(5.2) 



where —6 < i? n < 13 (we refer the reader to [8] for a more thorough discussion of both these results 
and more recent work on this topic). Returning to 0(n; 9) we now substitute the result of Lemma 
into (5.1) and (5.2) to obtain the global and asymptotic bounds. In Figure [9] we compare these 



5.1 



bounds to their numerically computed values. The relative sharpness of such estimates is once more 
observed. 



5.2 Polynomial sampling 

The primary concern of this paper has been reconstruction from Fourier (or Fourier-like) samples. 
However, in several circumstances, most notably the spectral approximation of PDEs with discon- 
tinuous solutions [23, 27 , the problem arises where a piecewise analytic function has been sampled 
in an orthogonal polynomial basis. As previously noted, this approximation will converge slowly 
(and suffer from a Gibbs-type phenomenon near discontinuities), hence it is necessary to seek a new 
approximant with faster convergence. Whilst a version of the Gegenbauer reconstruction method 
has been developed for this task HH] , the advantages of the method proposed in this paper (see 



Section 3.6 ) make it a potential alternative to this existing approach. Hence, the purpose of this 
section is to give a brief overview of this application. 

It is beyond the scope of this paper to develop this example of the reconstruction procedure 
in its full generality. Instead, we consider only the recovery of a piecewise analytic function / : 
— 1, 1] — > K from its first m Legendre polynomial coefficients fj = (/, j = 0, . . . , m — 1, where 

we 



ipj = (j + \ ) 2 Pj (x) is the j th normalised Legendre polynomial. Proceeding as in Section 



3.5 



J \J ~ 2 / J ' 

assume that / has jump discontinuities as — 1 < Xi < . . . < X\ < 1, and seek an approximation of 
the form 

l n r -l i 
fn,m(x) = X] X] a r,j4>r,j(x), U = ^Il r , 
r=0 j=0 r=0 

where 4>r,j{ x ) = -^4>j{^r{x)), K r (x) = - 1, c r = \{x r +i —x r ) and {0 O , ■ • ■ ,<t>n-i} is a system 
of polynomials on [—1, 1]. Since / is piecewise analytic, we expect exponential convergence of f n ,m 
to /, provided m is sufficiently large in comparison to n. 

Aside from determining how large m must be in comparison to n for recovery, the main question 
remaining is that of implementation, i.e. how to compute the entries of the matrix U. This requires 
evaluation of the integrals 

ipj(x)(f> r: k(x) dx, j — , . . . , in — 1 , k = 0, ...,n — 1. 

Whenever the reconstruction functions 4> r ,k arise from Gegenbauer polynomials, these calculations 
can be done iteratively. For the sake of brevity, we will not describe this computation in full 
generality. Instead, wc consider only the situation where the functions <f> r ^ arise from Legendre 
polynomials, in which case we are required to compute the integrals 

/ P,-(a:)Pfc(A r (a;))da;, j = 0, . . . , m - 1, k = 0, . . . ,n - 1, r = 0, . . . , I. 

J x r 

We have 

Lemma 5.2. Let 



[ Pj(x)P k (cx + d)dx, j, k = 0,1,2,..., (5.3) 

J a 
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where ca + d = — 1 and cb + d = 1 . Then 



"o,o = b - a, Uj. = - [P j+1 (x) - xPj(x)] h x _ a , j = 1,2,..., 

2 2 2d 

ua,k — -6o,k, uik = —5xt i S ok , k = 0,1,2,..., 

c 6c c 



and, for j > 2 and fc > 1, 



(2j-l)(fc + l) 



- u 



cj(2k+l) hk+1 " cj(2k + 1) " J 

Proof. Recall the recurrence relation 

jP.Or) = (2j - 1)^-1 (a:) - (j - ljPj-afr). J = 2, 3, . . . , 



(2j - l)fe 



: u 7-i ,fc-i 



<2j - 1) i - 1 

. «j-i,k — «^_2,fc. (5.4) 



c.7 



(5.5) 



for Legendre polynomials [TJ chpt 22]. Substituting this into (5.3) gives 



i - 1 

xPj-i(x)Pk{cx + d) dx ; — Uj-2,k- 



Letting x i— > cx + d in ( 5.5 ) and rearranging, we find that 
xP k (x 



cC&Tif^ + d) " ^ Pfc(CT + d) + c(2¥TT) Pfe - l(c * + 



The recurrence (5.4) now follows upon substituting this into the previous expression. 

Next consider ttj,o- Since Pq = 1, we have uo.o = b — a and u^o — / Pj{ x ) dx for j > > 1. Recall 
that the j th Legendre polynomial satisfies the Legendre differential equation [TJ chpt 22] 

[(l-z 2 )Pj(z)]' = -j(j + l)P j (x). 
Substituting for Pj in J Pj(x)dx and integrating gives 

The result now follows directly from the expression 

(1 - x 2 )P;(x) = (j + l)(xPj(x) - P j+ i(x)), j = 0, 1, 2, ... . 
To complete the proof, we consider uq^ and u\ tk - By the assumptions on a, h, c, d, we find that 

1 



Uo,k 



Pk(x) dx. 



Orthogonality now gives uq^ = 8o,k, as required. For u\ tk we have 



1 



Ul,k 



(y-d)P k (y)dy 



1 (2 



c V 3 



as required. 



□ 
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Figure 11: Left: the error log 10 \f(x) — f n ,m(x)\ for — 1 < x < 1 and m = 20,40,80, 160. Right: log errors 
l°Sio 11/ _ /n,m||oo (squares) and log 10 ||/ - fn,m\\ (circles) against m. 



n 


8 


16 


24 


32 


40 


48 


56 


64 


72 


80 


k(U*U) 


0.98 
19.97 


0.87 
4.17 


0.85 
3.57 


0.85 
3.57 


0.85 
3.50 


0.84 
3.43 


0.84 
3.43 


0.84 
3.41 


0.84 
3.38 


0.84 
3.38 



Table 5: The values C n , m and k(U*U) against n, where m = |n 2 . 



In Figure [IT] we consider the approximation of the function 



by the aforementioned method, using parameter values m = hn 2 , jiq — = \n and n\ — \n. As 



exhibited, we obtain 13 digits of accuracy using only m sa 120 Legendre coefficients of (5.6 1. Note 
that, although we have not shown it, the scaling m — O (n 2 ) appears to be sufficient for recovery. 
Numerical results, demonstrating this hypothesis, are given in Table [5} 



The function (5.6) was introduced in |28) to test the Gegenbauer reconstruction method when 
applied to this type of problem. As shown in Figure 11 we obtain a uniform error of roughly 10~ 8 



using only m = 40 coefficients, and when m — 120, the corresponding value is 10 -14 . In comparison, 
the Gegenbauer method gives errors of roughly 10~ 3 and 10 -7 for these values of m (see [28l Fig. 
3]), the latter being 10 7 times larger. 

Whilst this method appears to be a promising alternative, it should be mentioned that the recur- 
sive scheme to compute the entries of U requires O (m 2 ) operations. Since only O (run) operations 
are required to compute the approximation f n ,rm this is clearly less than ideal. Having said that, 
the Gegenbauer reconstruction method requires O (m 2 ) operations to compute each approximant, 
whereas with this scheme this higher cost is only incurred in a preprocessing step. 



6 Conclusions and future work 

We have presented a reconstruction procedure to recover a function using any collection of linearly 
independent vectors, given a finite number of its samples with respect to an arbitrary Riesz basis. 
This approach is both stable and quasi-optimal, provided the number of samples m is sufficiently 
large in comparison to the number of reconstruction vectors n. Moreover, this condition can be esti- 
mated numerically or, in certain circumstances, analytically. A prominent example of this approach 
is the reconstruction of a piecewise analytic function from its Fourier samples. Using a piecewise 
polynomial basis, this results in an approximation that converges root exponentially in terms of m, 
or exponentially in terms of n. 

There are many potential avenues to pursue in extending this work, as we now detail: 
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1. Piecewise polynomial reconstructions from polynomial samples. In the final section in the paper 
we detailed the recovery of a piecewise analytic function in a piecewise polynomial basis, given its 
Legendre polynomial expansion coefficients. Herein, an important open problem is verifying that 
the scaling m — O (n 2 ) is sufficient for reconstruction. Other challenges involve devising an iterative 
scheme for computing the entries of U valid for reconstructions in arbitrary Gegenbauer polynomials, 
and which involves only O (mn) operations, as opposed to O (to 2 ). Naturally, future work will also 
investigate the extension of this approach to reconstruction from arbitrary Gegenbauer polynomial 
expansion coefficients (as opposed to just Legendre polynomial expansion coefficients). 

2. Gegenbauer polynomial reconstructions from Fourier samples. As shown, the reconstruction 
procedure can be implemented with arbitrary Gegenbauer polynomials. However, unless Legendre 
polynomials are used, the reconstruction is not stable. This problem arises because Gegenbauer 
polynomials do not form a Riesz basis for the space L 2 (— 1, 1) unless A = |. However, Gegenbauer 
polynomial do form an orthogonal basis for the weighted space L 2 (— 1, 1), where u>(x) = (1 — x 2 ) A ~3 . 
Hence, it is natural to ask whether the reconstruction procedure can be adjusted to incorporate this 
additional structure, thereby yielding a stable method. It turns out that this can be done, with the 
first step being the derivation of an extended abstract framework along similar lines to Section [2] 
We are currently compiling results in this case, and will report the details in a future paper. 

3. Applications. Aside from the obvious applications in image and signal processing, the are a 
number of other potential uses of the procedure. First, it may be applicable to the spectral dis- 
cretisation of PDEs. Spectral methods are extremely efficient for solving problems with smooth 
solutions. However, for problems that develop discontinuities, e.g. hyperbolic conservation laws, a 
postprocessor is required to recover high accuracy [27] ■ The Gegenbauer reconstruction technique 
has recently been applied to such problems (see [231 HZ] and references therein) . Given the potential 



advantages of the method developed in this paper (see Section 3.6|, it is significant interest to apply 
this approach to these problems. Aside from high accuracy, a pertinent issue in the use of spec- 
tral approximations for nonsmooth problems is the question of stability |27j . Since the generalised 
reconstruction procedure is well-conditioned, there may also be a benefit in this regard. Outside 
of PDEs, the Gegenbauer reconstruction technique has also been extended to other types of series, 
including radial basis functions [39], Fourier-Bessel series [40] and spherical harmonics [22]. Future 
work will also consider generalisation of the method of this paper along these lines. 

4- Spline and wavelet-based reconstructions. Reconstructions in spline and wavelet bases are impor- 
tant in numerous applications. In [3], the authors gave a first insight into the application of such 
bases to the Fourier sampling problem. However, the theory is far from complete. Moreover, the 
use of more exotic objects, such a curvelets [E] and ridgelets [15], remains to be investigated. 

5. Recovery from pointwise samples. The discrete analogue of the Fourier coefficient recovery 
problem involves the reconstruction of a function from to equispaced samples in [—1, 1]. This problem 
has received more attention of late [13j [42] than the continuous case considered in this paper. In 
particular, the so-called mock-Chebyshev method [14] can be viewed as a discrete analogue of 
this approach. Whilst the mock-Chebyshev method is well understood, there remain a number of 
other reconstruction from pointwise samples problems of interest. In particular, with application 
to spectral collocation schemes based on Chebyshev or Legendre polynomials, the recovery of a 
piecewise analytic function from its values at Gauss or Gauss-Lobatto nodes. Future work will 
consider this problem. 
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